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in  the  combustion  chamber  at  any  Section.  The  gas  is  assumed  to  be  in 
chemical  equilibrium  and  is  coupled  to  the  spray  analysis  by  heat,  mass, 
and  momentum  transfer.  Combustion  gas  properties  are  evaluated  at  the 
local  temperature,  pressure,  and  O/F  ratio.  The  solution  of  the  problem 
requires  the  simultaneous  solution  of  six  ordinary  differential  equations. 

The  system  of  coupled  equations  are  integrated  and  solved  on  an  IBM  7040 
digital  computer.  The  analysis  is  applied  to  a  500  lbf.  nominal  thrust 
JP-5A,  liquid  oxygen  rocket  motor.  The  results  indicate  that  the  oxidizer 
evaporates  more  rapidly  than  the  fuel,  producing  a  combustion  gas  O/F  ratio 
during  the  initial  stages  of  evaporation  which  is  considerably  greater  than 
the  O/F  ratio  of  the  injected  propellants.  Therefore,  the  flame  temperature 
in  the  beginning  of  the  chamber,  and  the  subsequent  rate  of  fuel  evaporation 
are  less  than  those  predicted  by  a  model  which  assumes  proportional  rates 
of  oxidizer/fuel  evaporation.  , 

A  series  of  steady  state  and  wave  propagation  experiments  are  con¬ 
ducted  with  a  2-in.  -diam.  variable  length  rocket  motor  using  JP-5A  and 
liquid  oxygen.  A  simple  converging-diverging  nozzle  is  employed  to  yield 
high  chamber  exit  Mach  numbers  (0.  45).  The  steady  state  experiments  are 
primarily  concerned  with  the  measurement  of  axial  combustion  chamber 
pressure  gradients  in  order  to  correlate  the  spray  combustion  analysis  and 
provide  initial  conditions  for  wave  propagation  studies*,  The  data  obtained 
is  useful  in  establishing  design  criteria  ior  steady  state  operation  and 
optimization  of  the  rocket  chamber  geometry.  The  results  from  ten  different 
injector  configurations  indicate  that  the  slope  of  the  axial  pressure  gradient 
can  be  altered  by  changing  the  impingement  point  and/or  the  diameter  of  the 
fuel  orifice. 
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Wave  propagation  experiments  are  conducted  by  utilizing  the  com¬ 
bustion  chamber  as  a  driven  tube  and  mounting  a  diaphragm  and  high  pressure 
driver  tube  upstream  of  the  injector.  The  initial  results  indicate  that  the 
incident  wave  causes  drop  shattering,  and  prove  that  wavelet  generation, 
coalescence  and  wave  steepening  occur.  In  addition,  the  downstream  propa¬ 
gating  wave  is  shown  to  reflect  from  the  sonic  plane  causing  the  chamber- 
nozzle  cavity  to  behave  as  a  half-wave  resonator.  The  fundamental  mode 
of  longitudinal  oscillation  has  a  frequency  equal  to  that  predicted  from  simple 
acoustic  theory  with  mass  motion. 

Additional  experiments  are  conducted  to  determine  the  conditions  under 
which  input  disturbances  attenuate  or  amplify  in  liquid  propellant  rocket  motors. 
The  results  indicate  that  wave  steepening  and  pressure  amplification  are 
strongly  coupled  to  the  steady  state  gas  dynamic  flow  field  through  which  the 
wave  must  propagate.  Those  injector-chamber  configurations  which  result 
in  rapi"’  propellant  utilization  and  high  steady- state  axial  pressure  gradient 
slopes  tend  to  inhibit  wave  growth.  Less  rapid  conversion  to  gaseous  products, 
with  an  attendant  low  pressure  gradient,  provides  an  energy  and  mass  source 
to  drive  the  wave  and  amplify  the  base  pressure. 
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SECTION  I 

INTRODUCTION 

A.  Definition  of  the  Problem 

The  occurrence  of  high  frequency  osci.latory  combustion  is  of  major 
concern  in  the  development  of  rocket  motors.  Basically,  high  frequency 
instability  involves  finite  amplitude  pressure  wave  propagation  in  the  combus¬ 
tor  and  nozzle,  which  often  results  in  complete  destruction  of  the  engine.  A 
coupling  mechanism  between  the  wave  and  gas  dynamic  processes  supplies 
energy  to  sustain  the  oscillations.  The  instability  phenomena  has  been 
observed  in  liquid,  solid  and  pre-mixed  gaseous  rockets. 

The  mode  of  oscillation  is  primarily  dependent  upon  the  geometric 
configuration  of  the  chamber  nozzle  cavity.  Longitudinal  oscillations  are 
prevalent  in  thrust  chambers  with  large  length  to  diameter  ratios,  while 
transverse  modes  are  associated  with  large  diameter,  small  length  units. 

The  transverse  waves  manifest  themselves  as  radial,  tangential,  sloshing 
or  spinning  modes.  The  waveform  of  the  longitudinal  instability  can  be  either 
sinusoidal  or  shock  fronted  with  exponential  decay. 

Instabilities  may  be  initiated  as  a  result  of  a  spontaneous  action 
within  the  gas  dynamic  field  or  as  a  consequence  of  an  external  triggering 
action.  In  the  former  case,  small  disturbances  in  the  flow  are  amplified  by 
the  energy  addition  coupling  mechanism  with  the  combustion  process.  The 
triggering  action  type  of  instability  initiation  results  from  large  disturbances 
to  the  system,  such  as,  vibrations  from  the  engine  frame,  abrupt  changes  in 
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the  feed  syotem,  etc.  The  latter  mode  of  initiation  involves  nonlinear 
effects  at  the  very  beginning  of  the  instability.  Nonlinear  effects  are  pre¬ 
sent  with  any  finite  amplitude  wave  so  that  they  must  be  considered  even 
for  the  spontaneous  initiation  where  a  small  perturbation  grows  in  amplitude. 
The  nonlinearities  referred  to  here  manifest  themselves  ultimately  as  wave 
distortion  phenomena,  which  in  the  limit  cause  either  complete  attenuation 
or  shock  formation. 

Thuo  the  problem  of  combustion  instability  in  a  liquid  propellant 
rocket  engine  includes:  (1)  a  forcing  mechanism  which  initiates  oscillatory 
behavior;  (2)  an  energy  coupling  mechanism  which  sustains  these  oscilla¬ 
tions;  and  (3)  mechanisms  which  either  attenuate  or  amplify  the  oscillations. 
The  forcing  mechanisms  that  initiate  the  instability  are  generally  random 
in  nature  and  are  not  amenable  to  a  systematic  analysis.  The  solution  to 
the  instability  problem  therefore  resides  in  a  complete  understanding  of 
the  coupling  and  attenuation  mechanisms. 

Energy  coupling  between  the  propagating  wave  and  the  combustion 
process  is  dependent  upon  the  relative  magnitudes  of  characteristic  times 
associated  with  the  source  (combustion)  and  sink  (wave).  If  the  relaxation 
time  of  a  significant  transport  process  such  as  liquid  atomization,  vapori¬ 
zation  or  chemical  reaction  is  less  than  or  equal  to  the  wave  residence  time 
in  a  volume  element,  coupling  can  occur  with  resultant  unstable  operation. 
The  leading  edge  of  a  passing  wave  in  a  reacting  droplet  system  can  increase 
the  rate  of  evaporation,  which  couples  as  a  mass  source  to  the  trailing  edge 
of  the  passing  wave,  thus  producing  amplification.  It  can  be  further  seen 
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that,  for  a  long  relaxation  Lime,  the  mats  source  can  generate  wavelets 
that  coalesce  as  they  propagate  and  ultimately  overtake  the  initial  wave 
that  caused  the  disturbance,  again  causing  amplification.  In  addition,  as 
a  wave  propagates  in  a  gas,  it  deforms.  A  compression  wave  will  steepen 
in  a  decelerating  flow;  an  expansion  wave  v/ill  broaden  in  an  accelerating 
flow.  A  compression  wave  in  an  accelerating  flow  and  an  expansion  wave 
in  a  decelerating  flow  will  either  broaden  or  steepen,  depending  on  the  wave 
slope  and  the  velocity  gradient  of  the  gaseous  medium.  Therefore,  as  a 
wave  moves  in  a  rocket  chamber,  it  will  change  its  geometry,  which  then 
alters  its  residence  time  in  a  volume  element.  In  addition,  the  energy  in 
a  wave  is  a  function  of  its  velocity  and  waveform.  Thus,  the  wave  slope 
assumes  a  major  role  in  that  it  determines,  by  modifying  the  effective 
wavelength  and  amplitude,  the  nature  of  the  energy  or  mass  coupling  to  the 
propagating  wave  and  the  ultimate  stability  of  the  system. 

In  order  to  achieve  a  series  of  rational  engineering  principles  for 
designing  stable  rocket  engines  it  is  necessary  to  investigate  the  mechanisms 
that  cause  energy  coupling  and  wave  deformation.  This  thesis  is  directed 
to  that  end.  Section  II  presents  a  steady  state  aerothermochemical  analysis 
for  a  liquid  bipropellant  rocket  motor  with  small  contraction  ratio.  Both 
fuel  and  oxidizer  ballistics  are  analyzed  simultaneously.  The  solution  to 
this  problem  provides  a  means  for  synthesizing  injector- chamber-nozzle 
configurations  for  optimum  steady  state  operation,  as  well  as  providing  well 
defined  initial  and  boundary  conditions  for  wave  propagation  studies. 

Section  III  describes  the  experimental  facility  designed  for  this  investigation. 
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Section  IV  presents  the  results  of  an  experimental  investigation  of  steady 
state  behavior  in  a  liquid  bipropellant  rocket  motor.  In  addition,  longitu¬ 
dinal  wave  propagation  studies  were  conducted  in  order  to  define  the  para¬ 
meters  that  determine  whether  an  input  disturbance  will  attenuate  or  amplify. 
The  results  of  this  effort  are  also  presented  in  Section  IV. 

B.  Review  of  Previous  Spray  Combustion  Models 

A  review  of  the  various  analyses  proposed  for  describing  the  phe¬ 
nomenon  of  spray  combustion  in  liquid  propellant  rocket  motors  indicates 
that  all  are  built  upon  a  foundation  of  assumptions.  The  areas  of  concern 
included  in  these  assumptions  can  be  categorically  listed  as  follows; 

1.  Influence  of  chemical  kinetics 

2.  Vaporization  rate  limited  combustion 

3.  Injection  boundary  conditions 

4.  Relative  propellant  vaporization  rateB 

5.  The  effects  of  forced  convection 

6.  Chamber  pressure  variation 

7.  Droplet  shattering 

8.  Spray  distributions 

Several  investigations  have  treated  the  liquid  rocket  combustion 
problem  using  the  assumption  that  the  combustion  rate  is  controlled  by 
chemical  kinetics.  Miesse*  considered  a  single  chemical  reaction  of  first 
order  and  assumed  a  one-dimensional,  isothermal  system.  In  addition  it 
was  assumed  that  the  propellant  vaporization  could  be  described  as  a  linear 
regression  rate  of  the  droplet  surface  area  in  the  absence  of  convection. 
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The  implications  of  the  results  of  Miesse's  analysis  is  that  there  is  little 
or  no  difference  between  spray  combustion  or  premixed  gas  combustion. 

This  infers  that  chemical  kinetics  is  the  combustion  rate  limiting  process. 

2 

The  above  conclusions  conflict  with  those  of  Bittker  and  Brokaw 
who  calculated  maximum  chemical  space  heating  rates  in  combustion  pro¬ 
cesses.  The  results  represent  the  theoretical  maximum  possible  heat 
release  rate  to  be  expected  from  a  unit  volume  of  a  reacting  fuel-oxidant 
mixture.  Comparison  of  the  results  with  estimated  heat  release  rates  for 
the  same  propellants  in  experimental  rocket  engines  shows  that  the  latter 
are  several  orders  of  magnitude  less  than  the  calculated  maximum  chemical 
rate  for  a  given  propellant.  These  results  indicate  that  the  combustion 
process  in  conventional  liquid  rockets  is  much  less  rapid  than  it  would  be 
if  limited  by  chemical  kinetics.  Therefore,  the  great  bulk  of  the  work  in 
spray  combustion  is  not  influenced  by  chemical  kineticB,  but  rather  concludes 
that  the  physical  processes  of  liquid  atomization,  droplet  vaporization  and 

gas  phase  mixing  are  the  rate  controlling  parameters. 

3-8 

Most  spray  combustion  analysis,  are  based  upon  a  one -dimensional 
model,  i.  e.  ,  all  variables  are  functions  of  distance  from  the  injector  face. 
This  implies  the  assumption  that  gas  phase  mixing  is  instantaneous,  since 
any  realistic  description  of  the  phenomenon  requires  at  least  two  dimensions. 
The  processes  involved  in  liquid  atomization  are  indeed  not  uni-directional, 
so  that  droplet  vaporization  must  be  considered  as  the  dominant  process  in 


the  rate  controlling  mechanism.  Since  most  spray  combustion  models  are 
similar  in  the  manner  in  which  they  treat  the  vaporization  limited  combustion 
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gas  generation  rate,  it  being  assumed  that  the  combustion  effects  from 
single  droplets  are  additive,  a  detailed  discussion  of  each  model  would 
be  redundant.  Rather,  the  various  investigations  will  be  discussed  as  to 
the  manner  in  which  they  treat  the  various  assumptions  listed  previously. 

The  region  near  the  injector  face  has  long  been  a  problem  for  both 
the  rocket  designer  and  theoretician.  As  the  propellants  flow  from  the 
injector  in  the  form  of  ligaments  or  jets,  external  and  internal  forces  tend 
to  cause  atomization.  Large  transverse  gradients  of  oxidizer-fuel  ratio  are 
present  as  well  as  strong  recirculation  currents  of  combustion  gases.  The 
boundary  conditions  generally  employed  at  the  injector  face  are: 

1.  The  gas  velocity  is  zero  at  the  injector  face 

2.  The  liquid  propellant  sprays  enter  the  chamber  completely 
atomized  with  a  well  defined  droplet  size  distribution 

3.  The  initial  droplet  velocities  are  all  equal  to  the  liquid  stream 
injection  velocity 

Moot  investigators  actually  ignore  the  region  adjacent  to  the  injector 
face  by  applying  the  above  boundary  conditions  at  the  droplet  formation  or 
jet  break-up  distance  which  can  vary  from  1  to  4  inches  downstream, 
depending  on  the  type  of  injector  being  simulated.  A  notable  exception 
to  the  above  is  the  work  of  Lambiris,  (unpublished),  which  is  qualitatively 
described  in  Ref.  32.  According  to  this  model  the  combustion  chamber  is 
divided  into  two  regions:  (1)  the  injector  face  region  where  large  gradients 
exist  in  the  distribution  of  propellants,  and  (2)  the  remaining  portion  of  the 
chamber  and  nozzle  where  the  flow  of  liquids  and  gases  is  assumed  to  be 
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one-dimensional  and  uniform.  The  region  adjacent  to  the  injector  face 
is  divided  into  oxidizer  rich  zones,  fuel  rich  zones,  well  mixed  zones  and 
recirculation  zones,  v/here  no  liquid  is  present.  In  addition  an  attempt  it 
made  to  account  for  liquid  fragmentation  and  incomplete  atomization. 
Although  it  appears  that  this  is  a  three  space  dimension  transient  problem, 
the  author  reduces  it  to  a  quasi- steady  one-dimensional  problem  by  pre¬ 
scribing  all  of  the  zones  and  processes  as  functions  of  the  axial  distance. 

As  yet,  the  quantitative  description  of  some  of  the  above  processes  is  not 
known,  even  under  environmental  conditions  less  severe  than  those  encoun¬ 
tered  in  liquid  propellant  rocket  motors.  These  parameters,  therefore, 
must  be  studied  and  deduced  from  controlled  experiments  before  the  validity 

of  the  above  model  can  be  ascertained. 

9  10 

Some  investigations,  ’  ,  while  recognizing  that  differences  in 

relative  vaporization  rates  exist,  have  not  distinguished  between  fuel  and 

3  4  6  7 

oxidizer  droplets  in  the  analyses.  Still  other  models  .  ’  ’  ’  have  found 
it  convenient  to  assume  that  the  governing  equations  could  be  written  for 
the  fuel  spray  alone  and  that  the  oxidizer  and  fuel  vaporization  rates  ratio 
was  constant  and  equal  to  the  injected  mass  mixture  ratio.  This  assumption 
would  imply  a  relatively  constant  and  high  combustion  gas  temperature  since 
all  products  are  generated  as  the  result  of  adiabatic  chemical  combination 
at  a  near  stoichiometric  oxidizer-fuel  ratio.  In  addition,  this  assumption 
could  not  properly  model  a  thrust  chamber  operating  above  the  critical 
pressure  of  one  of  the  propellants. 
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The  vaporization  rates  of  both  propellant  species  must  be  con¬ 
sidered  in  estimating  the  overall  rate  of  combustion  product  formation. 


The  relative  rates  of  evaporation  of  fuel  and  oxidizer  depends  on  the  heat 
of  vaporization,  liquid  density,  thermal  capacity,  degree  of  atomization, 
the  critical  properties  of  the  propellants  and  the  operating  pressure  of  the 
thrust  chamber. 

7,  10 

Some  analytical  models,  ’  ,  are  based  on  the  assumption  of 

complete  entrainment  of  the  droplets  in  the  gas  stream.  In  the  absence 

of  relative  motion  between  the  gas  and  liquid  droplets,  forced  convective 

effects  augmenting  heat  and  mass  transfer  processes  are  neglected.  Another 

effect  of  forced  convection  on  a  droplet's  behavior  is  drag  with  resultant 

changes  in  droplet  ballistics.  Realistic  empirical  relations  for  Nusselt 

numbers  for  heat  and  mass  transfer  and  drag  coefficients  are  available 

3  4  6 

and  have  been  employed  in  several  investigations,  ’  ’  . 

As  combustion  proceeds  along  the  chamber,  the  stagnation  pressure 

will  decrease  and  the  gas  velocity  will  increase.  This  requires  that  the 

static  pressure  and  gas  density  decrease.  For  thrust  chambers  with  high 

contraction  ratios  the  decrease  in  static  pressure  is  negligible.  Many 

current  liquid  propellant  rocket  motors  have  small  contraction  ratios  such 

that  the  static  pressure  decrease  through  the  chamber  may  be  as  high  as 

20  percent.  The  decrease  in  static  pressure  causes  significant  increases 

in  vaporization  schedules  as  it  decreases  the  saturation  temperature  of  the 

propellant  and  increases  its  mass  diffusion  rate.  Most  spray  combustion 
3  6  7  9  10 

models  ’  ’  ’  ’  have  employed  an  assumption  of  constant  pressure 
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throughout  the  combustion  chamber.  Further  assuming  constant  molecu¬ 
lar  weight  for  the  gases,  coupled  with  the  previous  assumption  of  constant 
relative  vaporization  rates  and  hence  constant  gas  temperature,  implies 
that  the  combustion  gas  density  is  constant  throughout  the  chamber.  There¬ 
fore,  most  spray  combustion  models  have  analyzed  the  problem  of  a  liquid 

droplet  propagating  through  a  gas  with  constant  physical  properties.  Burstein 
4 

et  al.  have  accounted  for  the  pressure  variation  by  the  inclusion  of  an 
integrated  momentum  equation  and  coupling  of  the  droplet  ballistic  and  gas 
dynamic  equations. 

Aerodynamic  drag  forces  exerted  on  the  liquid  droplets  may  reach 
such  proportions  that  droplet  deformation  and  disintegration  will  eventu¬ 
ally  take  place.  Conditions  which  eventually  result  in  droplet  shattering 

11  12 

have  been  studied  experimentally  *  .  One  criterion  deduced  was  that 

breakup  can  be  expected  whenever  a  critical  value  of  the  Weber  number, 

(ratio  of  external  shear  forces  to  surface  tension  forces),  is  exceeded. 

After  the  critical  conditions  have  been  imposed  on  the  droplet,  shattering 

3 

may  be  delayed  and  non-uniform.  Lambiris,  et  al.  have  considered 
droplet  shattering  to  take  place  for  all  drops  over  50  microns  at  a  pre¬ 
scribed  distance  from  the  injector  face  based  upon  streak  photographs 

4 

obtained  from  a  transparent  rocket  motor.  Burstein  et  al.  have  included 
droplet  shattering  in  their  spray  combustion  model  by  replacing  the  shat¬ 
tered  droplets  by  an  equivalent  mass  of  smaller  droplets  whenever  a  pre¬ 
scribed  critical  Weber  number  was  exceeded. 

Spray  combustion  models  of  necessity,  must  take  into  account  the 
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number  of  droplets  that  constitute  the  spray  and  the  distribution  of  droplet 

sizes  among  that  number.  Extensive  reviews  of  the  statistical  description 
13  14 

of  sprays,  ’  ,  are  available.  The  log-probability  equation,  the  Rosin- 

Rammler  equation  and  the  Nukiyama- Tanasawa  equation  have  been  found 
to  fit  observed  distributions  reasonably  well.  Rather  than  employ  a  cumber¬ 
some  mathematical  expression  for  the  spray  distribution,  all  spray  combus¬ 
tion  models  utilize  either  a  single  mean  droplet  size  or  approximate  the 
spray  distribution  with  a  finite  number  of  droplet  size  groups. 

Consideration  of  the  analyses  and  assumptions  discussed  previously 
leads  to  the  following  conclusions  regarding  the  development  of  a  more 
realistic  model  for  spray  combustion  in  a  liquid  propellant  rocket  motor. 

1.  Spray  combustion  is  vaporization  rate  limited  and  the  effects 
of  chemical  kinetics  can  be  neglected. 

Z.  The  region  near  the  injector  face  is  not  amenable  to  incorporation 
in  a  one -dimensional  steady  state  combustion  model. 

3.  The  overall  spray  combustion  model  can  be  based  on  the  sum¬ 
mation  of  the  histories  of  a  large  number  of  individual  droplets. 

The  accuracy  of  the  final  results  is  more  dependent  on  the  cor¬ 
relation  between  the  assumed  droplet  sizes  and  the  actual  spray 
distribution  than  any  other  single  factor  in  the  analysis. 

4.  The  effects  of  forced  convection  on  heat  and  mass  transfer 
processes  must  not  be  neglected.  In  addition,  heat  transfer 
between  the  gas  and  liquid  must  be  carefully  considered  with 
particular  attention  focused  on  variations  in  droplet  temperature 


10 
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due  to  heating  up  periods  and  cooling  due  to  large  mass  dif¬ 
fusion  rates. 

5.  Variations  in  gas  pressure  and  other  properties  must  be 
included  in  the  analysis  of  small  contraction  ratio  motors. 

6.  Provisions  should  be  made  in  the  model  for  incorporating 
droplet  break-up  processes.  At  the  present  time  the  availa¬ 
bility  of  empirical  data  on  droplet  shattering  in  rocket  motor 
environments  is  sparse. 

7.  Vaporization  of  both  the  fuel  and  oxidizer  must  be  analyzed. 

The  implications  of  incorporating  this  item  in  a  spray  combustion 
model  are  worthy  of  discussion.  The  energy  added  to  the  gas 
stream  is  a  function  of  the  ratio  of  oxidizer  to  fuel  evaporated 
at  the  volume  element  being  considered.  The  local  gas  prop¬ 
erties,  (specific  heat,  enthalpy,  temperature,  molecular  weight), 
are  functions  of  the  ratio  of  oxidizer  to  fuel  evaporated  between 
the  injector  and  the  volume  element  being  considered.  These 
parameters  are  predicted  as  the  results  of  thermodynamic  equi¬ 
librium  calculations.  For  a  particular  propellant  combination 
wherein  one  specie  vaporizes  much  more  readily  than  the  other, 

f 

(either  due  to  high  volatility  or  operation  above  the.  critical  pres¬ 
sure),  the  gas  temperatures  at  the  upstream  end  of  the  chamber 
will  be  considerably  lower  than  those  calculated  for  combustion 
at  the  injected  O/F  ratio. 

A  model  of  bi-propellant  dpray  combustion  (considering  both  fuel 

11 
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and  oxidize r  evaporation)  accounting  for  all  of  the  assumptions  enumerated 
above  is  developed  in  Section  II. 

C.  Review  of  Previous  Experimental  Work 

This  section  is  concerned  with  a  review  of  experimental  investigations 
of  combustion  instability  in  liquid  propellant  rocket  motors  with  particular 
emphasis  on  diagnostic  studies  of  wave  propagation  phenomena.  Instability 
in  the  liquid  propellant  rocket  system  manifests  itself  as  uncontrolled  cyclic 
variations  of  pressure  in  the  feed  system  and  combustion  chamber  concur¬ 
rently  or  in  the  combustion  chamber  alone,  depending  upon  the  frequency 
of  oscillations.  High  order  frequencies,  greater  than  1000  cps  are  almost 
always  confined  to  the  combustion  chamber,  and  oscillations  having  fre¬ 
quencies  below  &00  cps  reflect  simultaneous  variations  in  the  feed  system. 

Initial  efforts  in  analyzing  low  frequency  instability  were  reported 

by  Gunder  and  Friant,  Subsequent  work  by  Summerfield,  Crocco, 17, 

18  19  EO 

Bar  re  re  and  Moutet,  ,  and  others,  ’  ,  have  shown  that  the  low  fre¬ 

quency  instabilities  may  be  related  to  a  difference  in  phase  between  the 
variations  of  the  chamber  pressure  and  the  burnt  flow  of  propellant  pro¬ 
duced  by  combustion.  This  phase  difference  results  from  the  time  lag 
between  injection  and  combustion  of  a  given  propellant  particle.  This  igni¬ 
tion  delay  was  at  first  considered  to  be  constant,  ^  ,  but  subsequently 

17 

was  related  to  the  fluctuations  of  the  pressure  in  the  chamber,  ,  The 
theoretical  analysis  predicted  stability  could  be  improved  by  increases  in 
the  length  of  the  feed  line,  the  propellant  velocity  in  the  feed  line,  the  pres¬ 
sure  drop  across  the  injector  and  the  ratio  of  chamber  volume  to  nozzle  area. 

12 
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In  an  exhaustive  experimental  program,  Barrere  and  Moutet,  cor¬ 
roborated  these  results  ;  d  found  that  the  frequency  of  oscillations  increase 
with  chamber  pressure,  decrease  with  characteristic  length,  (ratio  of 
chamber  volume  to  nozzle  area),  and  increase  with  injection  pressure 
drop.  In  general,  low  frequency  instabilities  are  ascribed  to  interactions 
between  the  processes  in  the  combustion  chamber  and  the  propellant  feed 
system.  It  is  sufficiently  understood  so  that  it  is  comparitively  easy  to 
avoid  or  cure. 

High  frequency  oscillatory  combustion  differs  from  the  low  frequency 
type  in  many  respects  besides  the  magnitude  of  the  observed  frequency  of 
combustion  chamber  pressure  fluctuations.  The  high  frequency  oscillations, 
or  "screaming",  are  characterized  by  high  amplitude  fluctuations  of  the 
chamber  pressure  only,  and  are  not  accompanied  by  pulsations  in  the  pro¬ 
pellant  flow.  In  addition,  screaming  results  in  a  significant  increase  in 
the  amount  of  heat  transferred  from  the  combustion  gases  to  the  walls  of 
the  chamber, whereas  low  frequency  instability  does  not  alter  the  heat  trans¬ 
fer  rates. 

Screaming  has  been  attributed  to  combustion-reinforced  pressure 
waves  passing  through  the  chamber  and  reflecting  from  the  chamber  sur¬ 
faces  to  trigger  succeeding  combustion  surges.  The  frequency  of  Buch 
waves  would  therefore  be  governed  by  the  velocity  of  wave  propagation  and 
the  geometry  of  the  chamber.  Although  the  equations  of  pure  acoustics  are 
insufficient  to  accurately  describe  the  wave  propagation  phenomena  in  an 
inhomogeneous,  accelerating,  chemically  reacting  flow,  the  oscillation 
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frequencies  may  be  expected  to  correspond  roughly  to  one  of  several 

simple  acoustical  modes  of  the  chamber. 

21 

Levine  and  Lawhead,  ,  showed  that  liquid  propellant  rocket 

motors  are  indeed  capable  of  sustaining  organ  pipe  type  oscillations, 

22 

Heidmann  and  Priem,  ,  made  pyrometry  studies  with  a  2  inch  diameter, 

27  inch  long  chamber  firing  liquid  oxygen  and  hydrocarbon  fuel.  They 

observed  the  two  lowest  modes  of  an  organ  pipe  oscillation  for  a  closed-cnd, 

open-end  system.  The  most  convincing  proof  of  the  existence  of  acoustic 

23 

type  waves  in  the  chamber  was  obtained  by  Berman  and  Logan,  ,  Berman 
24 

and  Cheney,  ,  and  the  staff  at  the  California  Institute  of  Technology's  Jet 

25 

Propulsion  Laboratory,  .  These  groups  have  published  excellent  streak 
photographs  of  the  longitudinal  waves  in  rocket  motors.  Their  experimental 
information  established  definitely  that  the  vibrations  of  the  combustion  gases 
occurring  during  oscillatory  combustion  correspond  to  the  propagation  of 
finite  pressure  disturbances  in  the  combustion  chamber.  Berman  and  his 
co-workers  found  that  the  pressure  waves,  which  appear  to  originate  near 
the  injector  in  the  combustion  zone,  are  initially  of  almost  perfect  sinu¬ 
soidal  shape.  After  passing  a  few  times  back  and  forth  through  the  chamber, 
the  waves  may  become  shocks  with  an  attendant  increase  in  peak  to  peak 
amplitude.  It  is  also  observed  that,  once  the  pressure  waves  travelling 
back  and  forth  in  the  chamber  assume  their  saw-toothed  shape,  the  vibrations 
of  the  combustion  gases  are  self  sustained  and  of  constant  amplitude  and 
frequency.  Thus  the  vibrations  must  be  maintained  by  a  feedback  mechanism 
caused  by  the  interaction  between  waves  and  reactions  in  the  chamber,  since 
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without  energy  addition,  waves  of  finite  amplitude  would  be  attenuated  by 
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viscous  and  thermal  dissipation.  Berman  and  Logan  attribute  the  oscil¬ 
lation  reinforcement  to  intermittent  ignition  of  the  accumulated  propellants 
near  the  injector  by  the  compressive  impulse  of  each  longitudinal  wave. 
They  also  hypothesize  that  strong  waves  ignite  droplets  which  burn  in  the 


wake  of  the  pressure  v.ave  as  it  propagates  through  the  chamber. 
25 


Ellis  et  al.  ,  also  found  that  the  organ-pipe  oscillations  were  of 
a  travellingwave  rather  than  a  standing  wave  character.  A  strong  pulse 
started  at  the  injector,  traveled  the  length  of  the  chamber,  was  reflected 
and  traveled  back  to  the  head  end.  During  the  cycle,  the  amplitude  of  the 
pulse  decreased  continuously.  The  weak  pulse,  upon  reaching  the  head  end, 
triggered  a  new- strong  pulse. 


26,  27 

Crocco  and  his  co-workers  have  related  the  stability  of  the 


rocket  engine  to  the  exponent  of  the  pressure  sensitivity  of  the  combustion 
reaction.  They  assume  that  combustion  time  is  inversely  proportional  to 
a  power  "n"  of  the  chamber  pressure.  The  value  of  "n"  necessary  to  drive 
an  oscillation  then  indicates  the  likelihood  of  that  oscillation;  the  higher  the 
value  of  "n",  the  less  likely  that  the  oscillation  can  set  up.  The  mathe¬ 
matical  model  chosen  by  Crocco  is  based  on  the  linearization  of  the  equations 
of  motion  for  small  perturbations.  Implicit  in  the  model  is  the  assumption 


that  the  Mach  number  in  the  chamber  i3  much  less  than  unity  .  The  experi- 


27 


mental  verification  of  the  theory,  ,  was  obtained  with  rocket  motors  that 


had  high  contraction  ratios  and  hence  low  Mach  numbers.  However,  the 

23-25,28 


existence  of  phenomena  involving  strong  waves  of  large  amplitude. 
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indicates  that  the  vibrations  of  the  chamber  gases  which  take  place  during 

oscillatory  combustion  are  non-linear  and,  therefore,  cannot  be  described 

accurately  by  means  of  linearized  equations. 

Recent  experimental  efforts  have  been  concerned  primarily  with 

the  gross  effects  of  variations  in  injector  configuration,  nozzle  geometry, 

propellant  mixture  ratio,  etc. ,  on  the  limits  of  high  frequency  combustion 
29-31 

instability,  -  Basic  investigations  on  the  interactions  of  controlled 

pressure  disturbances  with  well  defined  steady  state  rocket  motor  flow 
fields  have  been  ignored.  The  experimental  portion  of  the  present  investi¬ 
gation,  detailed  in  Section  IV,  is  concerned  with  first  determining  the 
pressure,  velocity  and  combusted  gas  distributions  in  a  liquid  bi-propellant 
rocket  motor  and  then  investigating  the  interaction  of  a  controlled  pressure 
wave  with  this  field. 
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SECTION  II 

STEADY  STATE  BIPROPELLANT  SPRAY  COMBUSTION  MODEL 
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A.  Description  of  the  Phenomena 

The  steady  state  transformation  of  a  pair  of  injected  liquid  pro¬ 
pellants  into  hot  combustion  gases  in  a  rocket  engine  is  the  primary  con¬ 
cern  of  this  section.  An  excellent  qualitative  discussion  of  the  various 
processes  involved  has  been  published  in  Ref.  32.  In  describing  the  over¬ 
all  transformation  it  is  best  to  separate  the  various  processes  involved  in 
order  to  determine  the  rate-controlling  factors. 

The  generation  of  hot  combustion  gases  begins  with  the  injection 
of  liquid  propellants  into  a  thrust  chamber.  Depending  on  the  type  of  injector 
employed,  e.  g.  ,  like  on  like  impinging,  unlike  impinging  or  showerhead, 
spray  fans  are  formed,  inside  which  the  propellant  streams  break  into  liga¬ 
ments,  and  eventually  into  small  droplets.  The  system  of  injected  streams, 
liquid  ligamentB  and  droplets  is  entirely  enveloped  by  hot  combustion  gases. 
Heat  is  transferred  from  the  combustion  gases  to  the  liquid  propellants, 
thereby  increasing  their  temperature  and  causing  them  to  vaporize.  The 
gases  also  exert  aerodynamic  forces  on  the  liquid  fragments  increasing  their 

rates  of  atomization,  vaporization,  and  altering  their  axial  velocities.  In 

33 

addition,  it  has  been  shown,  ,  that  the  gases  near  the  injector  are  in  a 
dissociated  state  and  are  reactive  with  either  fuel  or  oxidizer  vapors. 

The  nature  of  the  injection  process  results  in  large  spatial  variations 
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of  the  distribution  and  fragmentation  of  the  liquid  propellants.  There¬ 
fore,  gradients  in  gas  temperature,  composition  and  pressure  will  exist 
in  the  region  near  the  injector.  However,  slightly  downstream  of  the 
injector  face,  interactions  between  the  liquids  and  gases  will  tend  to  equalize 
the  gradients,  so  that  for  steady  state  conditions  a  unique  distribution  of 
propellant  weight  frac  tions,  droplet  size,  and  gas  properties  can  be  assumed. 
The  resulting  distribution  is  primarily  imposed  by  the  geometrical  shape 
of  the  injector  and  its  operating  conditions.  The  preceding  discussion 
implies  that  a  well  mixed,  one-dimensional  combustion  gas  and  propellant 
mixture  cannot  be  assumed  near  the  injector  face.  Rather,  some  jet  break¬ 
up  or  droplet  formation  distance  must  be  assumed,  and  spray  combustion 
calculations  commenced  at  that  point.  The  magnitude  of  the  jet  break-up 
distance  is  generally  between  1  and  4  inches,  depending  upon  the  type  of 
injector  employed. 

In  the  region  downstream  of  the  droplet  formation  distance,  the 
large  transverse  concentration  gradients  diminish  until  they  become  neg¬ 
ligible.  In  this  region  the  injected  liquid  jets  have  been  completely  atomized 
into  well  defined  droplet  distributions  so  that  the  liquid  oxidizer-fuel  ratio 
is  nearly  constant  over  the  entire  chamber  cross  section.  In  addition,  the 
evolved  gases  are  uniform  over  a  given  chamber  cross  sectional  area. 

The  volumetric  flow  rate  of  the  liquid  propellants  is  very  small 
compared  to  that  of  the  gases  so  that  droplet  collisions  and  interference 
can  be  neglected.  Heat  is  transferred  from  the  hot  gases  to  the  propellants 
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causing  them  to  vaporize,  interdiffuse,  and  subsequently  react  with  the 
surrounding  gases.  The  resulting  accelerating  gaseous  flow  imposes 
aerodynamic  forces  on  the  droplets  which  causes  deformation  of  their 
shape.  Eventually,  the  induced  shearing  forces  may  overcome  the  droplets 
cohesive  force  breaking  them  into  smaller  fragments.  This  phenomenon  is 
referred  to  as  droplet  shattering  or  secondary  atomization.  When  all  of  the 
droplets  have  been  completely  vaporized  the  transformation  process  is 
complete,  and  the  subsequent  flow  problem  is  one  of  gas  dynamics. 

As  a  result  of  the  preceding  remarks  it  appears  that  the  transformation 
of  liquid  propellants  into  hot  combustion  gases  is  rate  controlled  by  the 
heat  and  mass  transfer  processes  under  forced  convection  conditions.  For 
some  propellant  combinations,  notably  the  hypergolics,  liquid  phase  chemical 
reactions  can  become  a  controlling  factor  and  should  be  considered  where 
applicable. 

B.  Theoretical  Analysis 

The  model  of  spray  combustion  in  a  rocket  motor  is  best  divided 
into  two  parts  which  are  strongly  coupled.  The  first  part  concerns  itself 
with  the  liquid  propellant  droplets,  while  the  second  part  describes  the 
combustion  gas  dynamics.  Coupling  is  obtained  by  virtue  of  mass  trans¬ 
formation,  energy  interchange  and  drag  forces.  Therefore,  a  set  of  equa¬ 
tions  will  be  derived  that  describe  the  fuel  droplet  history  and  a  similar  set 
will  be  used  to  describe  the  oxidizer  droplet  history.  A  third  set  of  equations 
will  be  required  to  model  the  combustion  gas  flow  field.  The  entire  system 
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of  equations  will  then  be  solved  simultaneously  with  appropriate  boundary 
conditions. 

Liouid  Droplet  Model 

The  heating  and  vaporization  of  fuel  or  oxidizer  is  assumed  to 
begin  at  the  point  in  the  chamber  at  which  liquid  droplets  are  formed. 

At  this  location  the  injected  propellant  streams  are  represented  by  a  drop¬ 
let  size  distribution.  The  particular  distribution  employed  is  dependent 
upon  the  type  of  injector  and  empirical  data.  The  subsequent  evaporation  of 
the  oxidizer  and  fuel  sprays  is  then  represented  by  a  summation  of  the 
evaporation  of  a  finite  number  of  representative  droplets.  It  is  further 
assumed  that  mixing  and  reaction  rates  are  fast  and  that  reacted  products 
are  formed  as  soon  as  the  propellants  are  vaporized. 

A  schematic  model  of  fuel  or  oxidant  drops  vaporizing  in  a  rocket 
engine  is  shown  in  Fig.  2. 1.  The  liquid  droplet  is  shown  at  position  x  in 
the  chamber  corresponding  to  time  t  and  an  increment  later  corresponding 
to  x  +  Ax  and  t  +  At.  In  the  interval,  the  drop  velocity  changes  by  Av  and  the 
gas  velocity  b’’  Au.  While  the  droplet  travels  through  the  increment  heat 
is  transferred  to  its  surface  at  a  rate  q^  and  mass  leaves  the  surface  at 
a  rate  w.  Therefore,  the  droplet  mass  and  radius  change  by  an  amount 
Am  and  Ar  respectively.  In  addition,  the  droplet  temperature,  which  is 
assumed  uniform  throughout  the  drop,  changes  by  an  amount  AT  . 

Xj 

The  entire  vaporization  process  in  divided  into  sufficiently  small 
increments  so  that  steady- state  mass,  momentum,  and  heat  transfer  equations 
are  applicable  within  each  interval.  One  dimensional  steady- state  flow  is 
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assumed  for  both  the  dropiet  and  bulk  gas  models.  The  droplet  density, 
surface  tension,  viscosity,  specific  heat,  vapor  pressure  and  heat  of 
vaporization  are  taken  as  functions  of  temperature.  The  heat  of  reaction 
of  fuel  and  oxidizer  is  evaluated  at  the  ratio  of  oxidizer  to  fuel  evaporated 
within  the  interval  under  consideration.  Properties  of  the  bulk  gas,  i.  e. , 
molecular  weight,  specific  heat,  enthalpy,  and  viscosity  are  evaluated  as 
functions  of  temperature  and  the  ratio  of  oxidizer  to  fuel  evaporated  between 
the  injector  and  point  of  interest. 

Mass  Transfer 

36 

Following  the  correlation  of  Ranz  and  Marshall,  ,  thermal  diffusion 
is  neglected  and  mass  transfer  results  from  a  driving  force  in  the  direction 
of  film  diffusion.  The  following  relationship  is  obtained: 


,1/2 


2K  RTr 

Nu.  =  — =  2  +  0.  6  (Sc)1  /  J(Re)J 
w  D  M 
v  v 


(1) 


where  Nu.  =Nusselt  number  for  mass  transfer 
w 


K  ^mass  transfer  coefficient,  lbm/lbf-sec 
g 

R  =  Universal  gas  constant,  ft  Ibf/lbm-mole  °R 

—  ...  o_  Teas  +  Tdrop 

T  =  film  temperature,  R;  — “ — 

r  =drop  radius,  ft 

M  =  molecular  weight  of  drop  vapor,  lbm/lbm-mole 

,  2 , 

Dv  =  diffusivity,  ft  /sec 

Sc  =  Schmidt  number 

Re  =  Reynolds  number 
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The  diffusivity,  D  .,  can- be  obtained  from  empirical  correlations  , 


or  from  equations  given  in  Bird,  et  al. , 
is  used  in  this  work: 
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The  latter  method 


(p  o  )^3{T  T  )5/12 
vpcapcb'  1  ca  cb'  a  _  b 

D  =  - ; -  —  T 

v  /....’d/2  p  R 


(2) 


where  p^^  =  critical  pressure  of  droplet  substance,  lbf/ft 


p^k  =  average  critical  pressure  of  film  constituents,  lbf/ft 


T  =  critical  temperature  of  droplet  substance,  R 
ca 


T  ,  =  average  critical  temperature  of  film  constituents,  R 
cb 


p  =  static  gas  pressure,  lbf/ft 
.-6 


a  =  3.  2  10 
b  =  1.823 

M  =  molecular  weight  of  droplet,  lbm/lbm-mole 
a 

“b  =  average  molecular  weight  of  film  constituents,  lbm/lbm-mole 


and 


tr  = 


(T  T  ,  ) 
ca  cb 


1/2 


The  rate  of  mass  diffusion  from  the  droplet  is  then  calculated  from 

(3) 


w  =  K  A ,  p  a 
g  drv 


where  w  =  rate  of  mass  diffusion,  lbm/sec 


A,  =  drop  surface  area,  ft 
a 


p^  =  drop  vapor  pressure,  lbf/ft 
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a  is  an  engineering  correction  factor  which  converts  the  equimolal 


diffusion  coefficient,  K  ,  (which  considers  diffusion  from  the  gas  to  the 

g 


droplet  to  be  equal  to  diffusion  from  the  vaporizing  material  to  the  sur¬ 


roundings)  to  a  unidirectional  diffusion  coefficient  which  only  considers 


diffusion  from  the  droplet  to  the  environment. 


Heat  Transfer 


The  heat  transfer  to  the  drop  analysis  is  that  given  in  Refs.  4  and 


6.  The  total  heat  transferred  from  the  gases  q^,  goes  to  heat  the  liquid 


drop  q^ ,  vaporize  the  diffusing  vapor  q^,  and  heat  the  diffusing  vapor  "q^. 


The  heat  arriving  at  the  droplet  surface  q  ,  equals  the  sum  of  q  and  q  and 

V  Sj  K 


i6  given  by; 


\  =  hVVT*> z 


2  o 

where  h  =  convective  heat  transfer  coefficient,  BTU/ft  -sec-  F 


T^  =  gas  temperature,  R 


T  =  droplet  temperature,  R 

JU 


In  Eq.  (4),  Z  represents  the  ratio  of  the  heat  that  is  conducted  to  the 


droplet  surface  with  convection  and  mass  transfer,  to  the  heat  transfer  with 


pure  convection  and  no  mass  transfer. 


The  film  coefficient  is  determined  from  the  correlation  of  Ranz 


“V  *  ^  -■  tr- 

’ 1  ■■  )\  *  *  - 


t 


mm 


S 
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and  Marshall, 
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Ni^  =  ~  =  2+0.  6<Pr)1/3(Re)1/2 


(5) 


where  Nu^  =  Nusselt  number  for  heat  transfer 


=  mean  film  conductivity,  BTU/ft-sec-  F 


Pr  =  Prandtl  number,  vapor  film 


Re  ~  Reynolds  number, 


2rlu-v  p 


u  =  gas  velocity,  ft/ sec 


v  =  drop  velocity,  ft/ sec 

3 


,  p  =  gas  density,  lbm/ft‘ 

All  transport  properties,  dimensionless  numbers  and  density  for 
the  gas  are  evaluated  at  the  average  film  temperature,  T. 

Momentum  Transfer 


Aerodynamic  drag  is  the  mechanism  by  which  momentum  is  trans¬ 
ferred  between  the  droplets  and  gaseous  environment.  The  drag  forces  will 
either  decelerate  or  accelerate  the  droplets  as  the  velocity  of  the  drop 
approaches  that  of  the  surrounding  gases.  For  a  spherical  drop  immersed 


in  a  moving  fluid  the  drag  force  F^  is  given  by: 


m . 


d  dv 

F,  = - —  =  C_ 

d  g  dt  D  g 


where  the  drag  coefficient,  C^,  is  given  as. 


(u-v)  |  u-v| 
2 
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(6) 


CD  =  27  (Re) 


■0.  84 


(7) 
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Changes  in  Droplet  Properties 

Changes  in  droplet  properties,  (Ar,Av,AT  ),  can  be  calculated 

X/ 

from  the  set  of  preceding  equations  assuming  one  knows  the  instantaneous 
drop  radius,  drop  velocity,  drop  temperature  and  state  of  the  surrounding 
gaseous  environment.  The  change  in  drop  radius  with  time  is  determined 
from  the  mass  transfer  equations  along  with  the  continuity  equation  for  the 
liquid  drop; 


)  =  -w 


(8) 


or 


dr 

dt 


w 


P£Ad 


(9) 


The  change  in  drop  temperature  with  time  is  determined  from  an 
energy  balance  at  the  droplet  surface; 


where 


V9X 


and  is  the  heat  that  changes  the  drop  temperature.  qv  and  w  are  deter¬ 
mined  from  Eqs.  (4)  and  (3)  respectively,  while  X,  the  heat  of  vaporization 
is  a  function  of  the  drop  temperature.  By  assuming  that  the  droplet  tem¬ 
perature  is  uniform,  the  time  derivative  of  the  drop  temperature  is  related 
to  the  heat  transfer  as  follows ; 


q£  =  mCp£ 


or 


dt  me 

p£ 


(qv-w  X) 


(10) 
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Examination  of  Eq,  (10)  shows  that  if  the  heat  transfer  rate  q 

v 

is  greater  than  the  rate  ai  which  energy  is  carried  away  by  the  vaporizing 
liquid  w  X,  the  droplet  temperature  increases.  Referring  to  Eq.  (9)  it 
can  be  seen  that  it  is  possible  for  the  drop  radium  to  increase.  This 

generally  occurs  during  the  initial  stages  of  evaporation  when  q  is  large, 

dT .  v 

.  a 

w  is  small  and  — ^  is  positive.  As  the  rate  of  evaporation  increases,  . 

wX  eventually  exceeds  q^  and  the  drop  temperature  decreases.  This  in 
turn  causes  a  decrease  in  the  drop  vapor  pressure  which  subsequently 
decreases  the  evaporation  rate,  Eq.  (3).  Therefore,  a  balance  is  reached 
wherein  the  heat  transfer  rate  equals  the  heat  carried  away  by  the  vapor  and 
the  drop  temperature  remains  constant. 

The  change  in  drop  velocity  is  obtained  from  Eqs.  (6)  and  (7)  and 
is  given  by: 


dv  =  1  _£  ,lu-_yj(u-vj 

dt  8  D  p  r 


(11) 


The  previous  system  of  equations  are  used  for  both  fuel  drop  histories 
and  oxidizer  drop  histories  with  the  substitution  of  appropriate  functions  for 
the  physical  properties  of  the  droplet.  The  equations  are  solved  for  each 
of  the  n  sizes  of  fuel  drops  and  m  sizes  of  oxidizer  drops  that  are  assumed 
to  represent  the  injected  spray  distribution  of  each  propellant.  The  average 
critical  properties  of  the  film  constituents  referred  to  in  the  expression  for 
the  diffusivity,  D^,  Eq. (2),  are  taken  as  the  average  for  CO^  and  H^O  for 
the  JP-5A-liquid  oxygen  system  discussed  in  the  analysis.  The  dimensionless 
quantities,  Re,  Sc,  Pr  and  the  gas  density  are  evaluated  at  the  film 
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temperature  T,  defined  previously. 
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where  A  =  chamber  cross  sectional  area,  ft 
cp  =  rate  of  flow  of  gases,  lbm/sec 
The  integrated  continuity  equation  is  simply: 

(puA)  =  (cp) 


X 


(14) 


The  rate  of  gas  flow  cp  is  obtained  from  the  droplet  model  by  adding  the 


evaporation  of  all  of  the  fuel  and  oxidizer  drops,  as  follows; 


n+m  x  w  N 

vl  I  -5-*  * 

p=l  o  P 


(15) 


where  p=l--~----n  represents  fuel  drops, 

p  =  (n  +  1)  -  -  -  -  (n  +  m)  represents  oxidizer  drops, 


and  subscript  p  indicates  a  particular  size  (class)  fuel  or  oxidizer  droplet. 


th 


represents  the  number  of  drops  of  the  p —  class  passing  location  x  per 


second,  v  is  their  velocity  in  feet  per  second,  w  is  obtained  from  Eq.  (3). 
P  P 


The  steady  state  one -dimensional  momentum  equation  is 


i  d£  =  __u_  du  _i_  dee  - 

p  d x  g  dx  pAg  dx'  ' 


(16) 


where  F  is  the  force  per  pound  of  gas  due  to  drag  between  the  gas  and  the 
uroplets  and  the  last  term  on  the  right  hand  side  represents  the  force  due  to 
the  momentum  increase  of  the  gaseous  mass  generated  within  the  control 
volume.  The  drop  velocity  v  is  taken  as  the  average  of  the  velocities  of  the 
various  size  fuel  and  oxidizer  drops.  The  average  is  weighted  according  to 
the  amount  of  mass  evaporated  from  each  class  within  the  interval.  The  force 
F  is  given  by: 
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1 


n+m 

*f  =  y  r 

L  p  pA 
p=i 


(BFOR) 


(17) 


where  the  force  on  an  individual  drop  class  is: 

m  /dv  \ 

£  _ _ E_  ( _ E]  Nd 

P  PAg  c\dt/  f 

and  m^  is  the  mass  of  the  size  (class)  drop. 

Integrating  the  momentum  equation  yields; 


(18) 


^xUx+Ax  .  . 

Px+Ax  Px  g  Ux+Ax  Ux 


(BFOR) 


x+Ax 


(u-v) 


Ax  - 


x+Ax  /dcp 


A  ““  Ag  \  dx 

“c 


(19) 


For  the  purpose  of  the  above  integration  the  gas  density  is  assumed  to  be 
constant. 


The  steady  state  energy  equation  for  a  control  volume  is 

2\1 


jd_ 

dx 


*lh+T 


Q- W  +  h 


dco 

stdx 


(20) 


where  the  left  hand  side  of  Eq.  (20)  represents  the  change  in  the  energy 
carried  by  the  gas  across  the  control  surfaces  and 

Q  =  net  heat  transfer  to  the  gas  within  the  control  volume,  BTU/ft-sec 
W  =  work  done  by  the  gas  within  the  control  volume,  BTU/ft-sec 

h  =  stagnation  enthalpy  of  the  propellants  evaporated  within  the 
81 

control  volume,  BTU/lbin 

The  last  term  on  the  right  hand  side  of"  Eq.  (20)  represents  the 
energy  added  to  the  bulk  gas  flow  by  the  liquid  propellants  evaporated  within 
the  control  volume.  The  net  heat  transfer  to  the  gas  Q  is  related  to  the 


heat  transfer  to  the  drop  q  ,  defined  by  Eq.  (4). 
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of  the  combustion  gases,  which  appears  in  Eq.  (12)  is  also  a  function 
of  gas  temperature  and  the  O/F  ratio  of  the  combustion  gases.  The  stag¬ 
nation  enthalpy  of  the  evaporated  propellants,  h^is  dependent  upon  the 
drop  temperature,  drop  velocity  and  ratio  of  oxidizer  to  fuel  evaporated 
between  x  and  x+Ax.  The  method  of  calculating  the  stagnation  enthalpy  or 
energy  addition  from  the  evaporated  propellants  is  discussed  in  detail  in 
a  subsequent  paragraph  under  the  Method  of  Computation. 

The  system  of  equations  that  must  be  solved  consists  of  Eqs.  (9), 

(10)  and  (11)  written  in  finite  difference  form  together  with  Eqs.  (12),  (14), 
(19),  and  (23).  The  first  three  equations  must  of  course  be  solved  for  each 
of  the  classes  of  both  fuel  and  oxidizer  drops. 

Method  of  Computation 

The  general  scheme  for  solution  of  the  previous  system  of  equations 
is  to  calculate  fuel  and  oxidizer  droplet  histories  based  upon  an  assumed 
gas  dynamic  flow  field.  The  hydrodyna  r.ic  equations  of  the  combustion 
gases  are  then  solved  based  upon  results  from  the  drop  calculations. 
Iteration  between  the  two  systems  of  coupled  equations  is  performed  until 
convergence  is  obtained.  The  problem  has  been  programmed  for  solution 
on  an  IBM  7040. 

The  Eqs. ,  (9, 10,  11)  that  describe  the  history,  (T  ,  v,  r)  of  a  fuel  or 

Xj 

oxidizer  drop  are  ordinary  differential  equations  that  can  be  solved  if  initial 

\ 

drop  conditions  and  values  of  gas  temperature,  gas  pressure,  gas  velocity 
and  combustion  gas  O/F  ratio  are  known  at  all  points.  The  required  droplet 
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initial  conditions  are;  injected  mass  flow  rate  of  oxidizer,  injected  mass 
flow  rate  of  fuel,  number  of  dif'urent  size  drops  that  represent  the  fuel 
or  oxidizer  droplet  dietrioutions,  the  percentage  of  total  flow  of  fuel  or 
oxidizer  represented  by  each  of  the  drop  sizes,  droplet  formation  distances 
for  each  of  the  drop  sizes  and  initial  values  of  temperature,  radius,  and 
velocity  for  each  of  the  drop  sizes.  The  last  four  parameters  may  have 
different  values  for  each  of  the  drop  sizes  considered.  With  the  above 
conditions  specified. equations  (9,  10, 11)  are  solved  repeatedly  to  obtain 
T  (x),  r(x)  and  v(x)  of  a  single  drop  for  each  of  the  drop  sizes.  The  results 
are  then  weighted  according  to  the  mass  fraction  of  each  drop  size  and  then 
summed  over  all  drop  sizes  to  obtain  the  total  evaporation,  heat  transfer 
from  the  gas  and  drag  forces  between  the  liquid  and  gas.  Physical  prop¬ 
erties  required  to  perform  the  calculations  are  given  in  appropriate  sub¬ 
routines  containing  equations  or  tabular  data  with  interpolation  routines. 
Liquid  propellant  properties  are  evaluated  at  the  local  drop  temperature. 
Combustion  gas  properties  are  evaluated  at  the  local  gas  temperature  and 
local  combusted  gas  O/F  ratio. 

The  chamber  and  converging  nozzle  were  divided  into  128  differential 
segments  of  equal  length.  Each  of  these  segments  was  then  divided  into 
8  subdivisions  of  equal  length.  The  calculation  of  changes  in  droplet  pro¬ 
perties  is  then  performed  throughout  1024  increments.  However,  if  the 
droplet  temperature  changes  by  more  than  10°F  within  any  subdivision, the 
number  of  subdivisions  within  that  segment  is  doubled.  If  the  liquid 
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temperature  change  still  exceeds  10  F  the  number  of  subdivisions  is 
again  multiplied  by  two.  Reduction  of  the  sub- step  length  continues  until 
the  maximum  drop  temperature  change  criteria  is  met. 

Equations  (12, 14, 19,  23)  represent  the  system  of  equations  to  be 
solved  in  the  gas  dynamic  subroutine.  The  chamber-nozzle  configuration 
is  divided  into  the  same  number  of  major  segments  and  minor  subdivisions 
as  in  the  liquid  droplet  calculations.  The  chamber  cross  sectional  area 


as  a  function  of  axial  distance  from  the  injector  is  calculated  by  a  subroutine 

\ 

\ 


program.  The  terms  that  couple  the  liquid  drop  equations  to  the  gas  equa¬ 
tions,  i.  e. ,  cp(x),  BHET(x),  BFOR(x),  BFOV(x)  and  the  combusted  gas  O/F 
ratio  as  a  function  of  x  are  available  from  the  evaporation  subroutine.  A 
separate  calculation,  employing  the  four  gas  system  equations  is  performed 
in  each  minor  subdivision  from  x  to  x+Ax.  The  method  of  solution  consists 
of  assuming  a  value  for  p(x+Ax)  and  solving  for  u(x+Ax)  from  Eq.  (14),  p(x+Ax) 
from  Eq.  (19)  and  h(x+£x)  from  Eq.  (23).  The  enthalpy,  h,  determines  the 
temperature  T  by  use  of  a  function  subroutine.  The  density  p,  corresponding 
to  the  last  calculated  values  of  pressure  and  temperature  is  calculated  from 
the  equation  of  state.  The  iterative  procedure  is  continued  by  using  the  new 
density  to  solve  for  velocity,  pressure  and  temperature.  When  convergence 
within  the  minor  subdivision  has  been  achieved,  calculations  are  started  in 
the  next  subdivision.  Calculations  proceed  to  the  nozzle  throat  at  which 
point  the  pressure,  temperature  and  velocity  profiles,  p(x),  T(x)  and  u(x) 
are  compared  with  their  values  from  the  preceding  iteration.  If  convergence 
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has  not  been  attained,  the  liquid  droplet  calculations  are  performed  using  the 

last  set  of  gas  dynamic  properties.  This  in  turn  will  cause  the  gas  dynamic 

history  to  be  recalculated,  etc.  The  problem  is  considered  to  be  solved 

when  convergence  is  obtained  in  the  gas  dynamic  profiles. 

Properties  of  the  combustion  products  of  JP-5A  and  liquid  oxygen, 

(enthalpy,  molecular  weight  and  specific  heat),  were  obtained  from  chemical 
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equilibrium  calculations  by  Barber  and  Hersch.  '  These  calculations  tab¬ 
ulate  the  above  properties  as  functions  of  temperature  and  O/F  ratio  at 
different  pressure  levels.  In  the  present  spray  combustion  analysis  the 
enthalpy,  molecular  weight  and  specific  heat  are  evaluated  at  the  local 


temperature  and  combusted  gas  O/F  ratio  by  means  of  a  double  interpolation 
subroutine.  Equilibrium  temperature  can  similarly  be  obtained  by  means 


of  a  double  interpolation  of  enthalpy  and  O/F  ratio. 

The  enthalpy  of  the  gaseous  combustion  products  as  given  in  the 
equilibrium  calculations  is  based  on  datum  enthalpies  of  -653BTU/lbm  at 
536°R  for  the  liquid  fuel  and -173  BTU/lbm  at  162 °R  for  the  liquid  oxygen. 
Accordingly,  the  term  in  the  energy  equation  (23)  which  represents  the  energy 
addition  to  the  gas  stream  by  the  evaporated  propellants  takes  these  datum 


states  in  account. 


h  42  =  42 

st  dx  dx 


(el5?7t]  t'635  '  173<ELOFi]  +  BHES 


where  ELOF  represents  the  ratio  of  liquid  oxidizer  to  liquid  fuel  evaporated 
within  the  minor  subdivision  under  consideration.  The  term  BHES  accounts 
for  the  fact  that  the  propellants  are  not  at  the  datum  state  prior  to  reaction 


but  rather  are  in  a  vapor  state  at  an  elevated  temperature  moving  with 


a  finite  velocity. 
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where  p  =  1 - n  represents  fuel  drops 

p  =  n+l - n  +  m  represents  oxidizer  drops 

Tj  =  datum  temperature,  °P. 

T  =  temperature  of.  drop,  °R 
Ju 

X  =  heat  of  vaporization,  BTU/lbm 

Am  =  reduction  in  droplet  mass  within  the  subdivision  A  t 
P 

N  =  number  of  drops  per  second  of  the  p  class. 


(26) 


The  calculation  within  the  brackets  is  performed  for  each  drop  size  indi¬ 
vidually  and  the  summation  is  then  carried  out  over  all  the  classes  of  fuel 
and  oxidizer  drops. 

The  integral  form  of  the  gas  dynamic  conservation  equations  written 
in  finite  difference  form  and  the  equation  of  state  contain  only  two  indepen¬ 
dent  variables;  pressure  and  temperature.  Combustion  gas  flow  rate  and 
O/F  ratio  are  specified  by  the  results  of  the  droplet  calculations.  Gas  density 
is  calculated  from  the  equation  of  state  as  a  function  of  pressure,  tem¬ 
perature  and  molecular  weight,  which  in  turn  is  a  fvinction  of  temperature 
and  O/F  ratio.  Gas  velocity  is  determined  from  the  algebraic  form  of  the 
continuity  equation  once  the  gas  density  has  been  calculated.  The  boundary 
conditions  on  the  gas  dynamic  equations  are  therefore  given  as: 
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at  x  =  x.  p  =  p 
1  o 

T  =  T 

o 

where  x^  is  the  point  at  which  gas  dynamic  calculations  commence  ;  the 

break-up  point.  If  in  addition,  the  existence  of  a  sonic  plane  is  specified, 

i.  e. ,  Mach  number  equals  one  at  x  =  x  ,  the  problem  becomes  overly 

constrained  since  the  injected  mass  flow  rate  and  cross-sectional  area  have 

previously  been  specified.  Therefore,  the  constraint  on  pressure  at  x.  is 

removed  and  the  problem  is  treated  as  a  split  boundary  value  problem.  An 

initial  value  of  p  is  chosen  to  start  the  gas  dynamic  calculations.  If  the 
o 

Mach  number  becomes  unity  before  the  throat  or  if  M  is  less  than  one  at  the 

throat,the  value  of  p^  is  changed  and  the  gas  dynamic  calculations  are  repeated. 

The  true  boundary  conditions  for  the  problem  considered  here  are  therefore, 

x  =  x.  T  =  T 
\  o 

x  =  x  M  =  1 
L 

The  above  boundary  conditions  are  applicable  to  systems  wherein  a 
chamber-nozzle  configuration  is  specified  and  the  problem  is  basically  one 
of  analysis  of  spray  combustion.  The  pressure-temperature  boundary  con¬ 
dition  is  more  appropriate  for  synthesis  of  a  system  where  the  chamber  pres¬ 
sure,  chamber  diameter  and  injected  mass  flow  rate  are  specified  and  the 
optimum  chamber  length  and  nozzle  area  variation  are  to  be  determined. 

C.  Results 

The  previous  analysis  and  its  attendant  computer  program,  described 
in  the  Appendix,  are  applied  to  several  sets  of  boundary  conditions.  For  all 
cases  tested,  the  thrust  chamber  geometry  is  similar  to  that  used  in  the 
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experimental  program.  The  chamber  is  2  in.  in  diameter  by  1.2  ft.  long. 

The  nozzle  converges  in  0.  1  ft,  to  a  throat  diameter  of  1.  64  in.  The  injec¬ 
tion  rate  of  fuel  and  oxidizer  is  .  63  lbm/sec  and  1.  701  lbm/oec  respectively. 

In  order  to  test  the  computer  program  and  compare  it  with  previous 

results,  it  is  initially  modified  to  simulate  the  "proportional  evaporation" 

4 

model  of  Bur  stein,  et.  al.  This  is  accomplished  by  deleting  the  oxidizer 

calculations  and  determining  the  combustion  gas  flow  by  replacing  Eq.  15  with; 

x  w  N  \  / 
j  dxj  fl  +  (O/F)  i 


where  p  =  1. . . .  n  represents  the  fuel  drops 

(O/F)  inj  =  ratio  of  injected  oxidizer  to  injected  fuel. 

As  a  result  of  Eq.  (27)  the  O/F  ratio  of  the  combustion  gas  is  constant 
throughout  the  motor.  The  combustion  gas  properties  data  of  Reference  4  are 
substituted  for  the  chemical  equilibrium  data  to  be  used  subsequently  with 
the  two-component  model.  This  involves  using  a  frozen  composition  gas 
specific  heat  given  as  a  function  of  temperature  only,  and  substituting  the 
product  of  specific  heat  and  temperature  for  enthalpy  in  the  energy  equation. 
The  energy  contribution  of  the  evaporated  propellants,  h^,  is  modified  to  be 
consistent  with  -V  change  in  datum  temperature. 

The  fuel  spray  is  approximated  by  four  different  drop  size  groups 
which  are  formed  at  0.  086  feet  from  the  injector.  The  boundary  condtions 
on  the  spray  are; 
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radius  (ft. ) 

concentration 

temperature  (°R) 

velocity(fps) 

-4 

1(10  -) 

.  30 

540 

142 

2(10 "^) 

.  35 

540 

142 

-4 

3(10  -) 

.  30 

540 

142 

4(10‘4) 

.05 

540 

142 

The  boundary  conditions  on  the  gas  system  are; 

x  =  .086  ft.  T  =  5000  °R 

x  =  1.  3  ft.  M  =  1. 0 

The  results  of  the  calculations  are  shown  in  Figs.  2.  2,  2.  3  and  2.4.  The 
drop  radius  initially  increases  due  to  a  decrease  in  density  as  the  drop 
temperature  rises.  The  distance  required  to  completely  vaporize  a  drop 
is  inversely  proportional  to  the  initial  drop  radius.  The  liquid  temperature 
rises  rapidly  until  an  equilibrium  condition  is  achieved,  at  which  point  the 
evaporation  rate  becomes  a  maximum  and  the  drop  temperature  remains 
constant.  In  general,  the  results  are  in  complete  agreement  with  those  of 
Bur  stein,  et.  al.  ,  wherein  it  is  stated  that  the  phenomena  is  one  of  slowly 
evaporating  drops  flowing  down  a  field  of  gradually  increasing  velocity  until 
a  cooperative  effect  between  the  drops  and  flow  field  causes  acceleration  of 
evaporation. 

The  computer  program  is  next  applied  to  a  two  component  system  as 
described  in  Part  B  of  this  section.  The  boundary  conditions  on  the  fuel  are 
as  given  previously.  The  liquid  oxygen  boundary  conditions  are: 
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radiue(ft.  )  concentration  temperature(°R)  velocity(fps) 


1(10’’) 

.  10 

160 

100 

V" 

1 

O 
•— < 

.  20 

160 

100 

3(10'4) 

.40 

160 

100 

»— » 
o 

t 

.30 

160 

100 

The  gas  system  boundary  conditions  are: 

x  =  .  086  ft.  T  =  3000°R 

x  =  1. 3  M  =  1 

The  results  are  shown  in  Figs.  2.  5,  2.  6,  2.  7  and  2.  8.  The  more  rapid 
heat  up  and  evaporation  of  the  liquid  as  compared  to  the  results  of  Burstein, 
et.  al.  is  not  expected.  However,  the  phenomena  is  readily  accounted  for 
when  one  considers  the  variation  in  combustion  gas  property  data  used  for  the 
two  analyses.  The  frozen  composition  gas  specific  heat  of  Burstein  is  only 
25%  of  the  equilibrium  composition  specific  heat  used  in  the  present  analysis. 
Since  the  heat  transfer  coefficient  to  the  drops  is  directly  proportional  to  the 
specific  heat,  the  bi-propellant  spray  analysis  with  equilibrium  gas  composition 
results  in  a  more  rapid  heat-up  and  vaporization  process.  In  order  to  compare 
a  proportional  evaporation  model  with  a  bi-propellant  evaporation  model,  the 
computer  program  is  again  modified  to  delete  the  oxidizer  calculations. 
However,  the  combustion  gas  data  is  now  taken  for  equilibrium  composition. 

The  results  of  these  calculations  are  shown  in  Figs.  2.  5  and  2.  7.  Only  the 
results  for  the  largest  size  drop  considered  are  presented.  The  figures 
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clearly  indicate  that  a  bi- propellant  evaporation  model  results  in  longer 
heat-up  periods  and  vaporization  distances.  Moreover,  comparison  of 
the  results  in  Figs.  2.  2  -  2.  4  with  the  proportional  evaporation  results  in 
Figs.  2.  5  and  2.  7  indicate  the  important  effect  of  gas  composition  data. 

The  boundary  conditions  for  the  spray  drop  size  distribution  used 
in  the  previous  examples  were  arbitrary  selections.  In  order  to  simulate  the 
conditions  encountered  in  the  experimental  engine  described  in  Section  HI 
and  Section  IV  of  this  report,  drop  size  distributions  are  determined  fol¬ 
lowing  the  method  of  Priem,  (6).  The  following  spray  conditions  are  obtained: 
radius  (ft)  concentration 

1. 5{10'4}  .30 

4(1 0  ~  4)  .40 

9(10-4)  .30 

This  distribution  is  assumed  to  approximate  both  the  fuel  and  oxidizer  spray. 
In  iddition,  the  fuel  drops  are  all  assumed  to  have  an  initial  velocity  of 
142  ft/ sec  and  an  initial  temperature  of  540°R.  The  oxidizer  initial  tem¬ 
perature  is  160°R  and  the  injection  velocity  is  100  ft/ Bee.  The  boundary 
conditions  on  the  gas  system  are: 

x  =  .086  ft.  T  =  3000°R 
x  =  1.3  ft.  M  =  1 

The  results  of  the  calculations  are  shown  in  Figs.  2.  9,  2. 10  and  2.  11.  The 
effect  of  using  larger  drop  sizes  to  approximate  the  spray  distribution  is 
clearly  shown  by  comparison  with  Figs,  2.5-2.  8.  The  liquid  oxygen 
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evaporation  is  not  retarded  as  much  as  the  fuel  evaporation,  because  of  its 
very  high  vapor  pressure.  As  a  result  the  O/F  ratio  of  the  combustion  gas 
during  the  initial  stages  of  the  chamber  is  increased  from  approximately  10 
in  the  previous  case  to  30  in  the  present  example.  This  causes  the  gas 
temperature  to  decrease  initially,  which  in  turn  reduces  the  heat  transfer 
rate  to  the  drops. 

Figure  2.  9  shows  the  variation  in  drop  temperature  with  distance 
for  the  various  fuel  and  oxidizer  drop  sizes.  The  effect  of  the  increased 
gas  O/F  ratio  and  reduced  gas  temperature  is  evidenced  by  the  long  heat-up 
periods  required  for  the  larger  fuel  drops.  Figures  2.  10  and  2.  11  show  the 
variation  in  fuel  and  oxidizer  evaporation.  The  results  of  the  gas  system 
calculations  are  shown  in  Fig.  2.  12.  In  addition,  a  curve  of  pressure  vs. 
distance  is  presented  for  a  proportional  evaporation  model,  using  equilibrium 
gas  composition  and  the  present  spray  distribution.  The  proportional  evapora¬ 
tion  model,  by  definition,  yields  a  constant  gas  O/F  ratio  of  2.  7  and  hence 
a  gas  temperature  profile  similar  to  that  shown  in  Fig.  2.4.  Therefore,  the 
subsequent  rate  of  gas  evolution  and  acceleration  is  move  rapid  causing  the 
pressure  to  decrease  closer  to  the  injector  face. 

A  series  of  experimental  data  points,  which  fall  between  the  two 
theoretical  curves  is  also  presented  in  Fig.  2.  12.  It  is  expected  that  experi¬ 
mental  pressure  points  fall  below  a  theoretical  curve  as  a  result  of  heat 
transfer  to  the  motor  walls  and  frictional  effects  which  are  neglected  in  the 
analysis.  It  appears  that  in  terms  of  the  present  experimental  data,  the  bi¬ 
propellant  evaporation  model  predicts  an  evaporation  rate  that  is  too  low, 


while  the  proportional  evaporation  model  yields  a  mass  liberation  profile 

that  is  too  high.  Referring  to  Fig.  2.  8,  the  pressure  profile  obtained  with 

-4  -4 

a  drop  size  distribution  ranging  from  1(10  )  to  4(10  )  ft  is  also  below 

the  experimental  data  points.  Note  also,  that  the  rate  of  gas  evolution  with 
the  smaller  aize  drop  distribution,  Fig.  2.  8,  is  rapid  enough  to  prevent  the 
decrease  in  gas  temperature  that  is  obtained  with  the  larger  drop  size  distri¬ 
bution,  Fig.  2.  12. 

The  previous  results  indicate  the  importance  of  the  boundary  conditions 
in  terms  of  proper  modeling  of  a  rocket  combustor.  It  is  reasonable  to 
assume  that  a  realistic  fuel  drop  size  distribution  exists  that  would  yield 
a  pressure  curve  in  complete  agreement  with  the  data  points.  It  has 
been  stated  previously  that  the  drop  size  distribution  of  Fig.  2.  8  was 
chasen  arbitrarily  while  the  logarithmiconormal  distribution  of  Fig.  2.  12 
was  determined  by  a  method  suggested  by  Priem  (6).  This  latter  method 
uses  the  results  of  cold-flow  tests  of  an  injector  to  determine  a  preliminary 
mass -median  drop  radius  in  terms  of  a  jet  diameter  and  a  velocity  difference 
between  gas  and  liquid.  The  preliminary  mass  median  drop  radius  is  then 
modified  to  account  for  injection  into  a  hot  combustor.  The  method  of 
modification  involves  a  vaporization  calculation  as  proposed  by  Priem.  The 
assumptions  required  to  determine  a  mass  median  drop  radius  for  combustor 
calculations  using  Priems  analysis  are  not  minor,  and  as  a  result,  large 
deviations  can  occur  between  experimental  and  theoretical  results. 

The  boundary  condition  for  gas  temperature  was  selected  after  several 
trial  computer  solutions.  The  value  of  JOOO°R  corresponds  to  the  average 
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adiabatic  flame  temperature  for  the  0/F  ratio  of  the  combustion  gases  in 
the  first  interval  of  calculation.  While  this  v?lue  is  suspect,  it  has  been 
shown  that  when  a  small  size  drop  distribution  is  used,  Fig.  2.  8,  the  gas 
temperature  does  not  diminish,  indicating  that  this  boundary  condition  plays 
a  secondary  role.  However,  any  future  modification  of  the  bi-propellant 
spray  combustion  model  should  include  a  method  to  vary  the  gas  temperature 
at  the  boundary  in  accordance  with  the  O/F  ratio  calculated  from  the  previous 
iteration.  In  addition,  the  effects  of  recirculation,  which  would  tend  to 
increase  the  gas  temperature  in  the  injector  region  should  be  examined  for 
inclusion  in  the  analysis. 

A  more  sophisticated  analysis  of  droplet  ballistics  might  include  heat 
conduction  through  the  droplet,  with  a  subsequent  radial  temperature  distri¬ 
bution.  This  would  cause  the  vaporization  rate  to  be  increased,  since  the 
surface  temperature  and  hence  the  vapor  pressure  would  rise  more  rapidly. 
The  phenomena  would  be  significantly  more  predominant  for  the  fuel  drops, 
resulting  in  a  lower  gas  O/F  ratio  with  a  subsequent  increase  in  gas  tempera¬ 
ture.  Therefore  the  rate  of  gas  evolution  would  be  increased  leading  to  a 
decrease  in  the  gas  pressure. 
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SECTION  m 

EXPERIMENTAL  FACILITY 

The  experimental  facility  has  been  designed  with  a  view  towards 
obtaining  both  steady  state  and  wave  propagation  data  in  a  liquid  propellant 
rocket  system.  The  steady  state  data  is  acquired  and  used  for  verification 
of  theoretical  spray  combustion  models  and  as  boundary  and  initial  conditions 
for  theoretical  and  experimental  wave  propagation  studies.  Information 
obtained  from  wave  propagation  programs  is  used  to  determine  the  para¬ 
meters  which  influence  energy  coupling  to  the  wave,  wave  coalescence, 
and  the  effect  of  velocity  and  pressure  gradients  on  wave  attenuation. 

The  test  site  complex  is  composed  of  a  control  building,  an  instru¬ 
mentation  building  and  a  test  cell  building.  The  buildings  are  connected 
via  underground  conduits  which  carry  hydraulic  lines  and  electrical  cable. 

The  test  cell  building  houses  four  individual  test  cells  which  contain  a 
500  lbf  thrust  stand;  a  4000  lbf  thrust  stand;  a  shock  tube  for  instrumentation 
calibration  and  a  hydraulic  bench  for  injector,  valve  and  flowmeter  calibra¬ 
tion.  One  wall  of  the  test  building  is  composed  of  sliding  doors  which  are 
rolled  back  to  afford  visual  observation  from  the  control  building.  The 
control  building  houses  the  test  console  from  which  all  command  signals 
are  initiated.  Visual  observation  of  the  actual  firing  is  accomplished  through 
a  blastproof  window.  The  instrumentation  building  contains  all  of  the  re¬ 
cording  and  peripheral  equipment  necessary  for  data  acquisition. 

The  rocket  engine  system  used  in  this  investigation  is  of  500  lbf 
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nominal  thrust.  The  propellants  are  JP-5A  and  liquid  oxygen.  The  rocket 
motor  itself  consists  of  a  stainless  steel  injector,  brass  or  bearing  bronze 
combustion  chamber  and  a  stainless  steel  nozzle  in  separable  units,  as 
shown  in  Fig.  3.  1.  The  combustion  chamber  is  made  from  one  or  more 
standardized  sections  depending  upon  the  desired  area  schedule  and  length. 
The  three  basic  components,  i.  e.  ,  injector,  chamber  and  nozzle  are  held 
together  by  a  mechanical  clamping  arrangement  with  teflon  "o- rings"  used 
to  seal  the  joints.  The  motor  is  uncooled,  and  uses  the  relatively  large 
mass  of  metal  to  absorb  the  heat  transmitted  from  the  gases. 

The  injector  head  assembly  (Fig.  3.  2)  consists  of  an  injector  plate 
and  a  two  part  manifold.  The  plate  threads  into  the  rear  manifold  section 
which  is  then  bolted  to  the  front  manifold  section.  Teflon  "o- rings"  prevent 
leakage  of  fuel  and  oxidizer.  The  assembly  was  designed  in  a  modular 
fashion  in  order  to  facilitate  changing  injector  plate  configurations.  An 
injector  plate  blank  is  shown  in  Fig.  3.  3.  A  detailed  discussion  of  the 
injector  plate  configurations  used  in  the  investigation  is  presented  in 
Section  IV. 

The  design  of  the  nozzle  was  influenced  by  the  objectives  of  the 
program  at  the  Propulsion  Research  Laboratory.  Among  the  objectives 
is  an  investigation  of  the  nonlinear  aspects  of  wave  propagation  and  com¬ 
bustion  instability  in  liquid  propellant  rocket  motors.  In  the  theoretical 

analysis  of  combustion  instability,  the  effect  of  the  Mach  number  of  the 

27  n 

gases  is  highly  significant,  as  indicated  in  .  Briefly  stated,  if  (1-Mj 
is  approximately  unity,  the  mathematical  analysis  can  be  greatly  simplified 
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by  means  o£  linearizing  assumptions.  However,  the  actual  problem  is 
23  24  28 

nonlinear,  ’  '  ,  and  in  practice  the  Mach  numbers  at  the  nozzle 

entrance  are  well  above  those  consistent  with  linearized  analyses.  As 
a  result  of  the  preceding  remarks  the  nozzle  contraction  ratio  was  designed 
for  a  nozzle  entrance  Mach  number  of  0.45,  based  on  isentropic  flow  and 
constant  isentropic  exponent. 

The  verification  of  the  nozzle  entrance  Mach  number  required  a 
measurement  of  the  total  pressure  at  that  point.  This  entailed  designing 
a  probe  that  can  withstand  extremely  high  temperatures  and  a  reactive 
atmosphere.  A  water-film  cooled  graphite  probe  was  developed  and 
employed  successfully  (Fig.  3.4).  The  cooling  water  is  sprayed  through 
the  probe  body  into  the  gas  stream  at  a  pressure  slightly  in  excess  of  the 
gas  pressure.  Thus,  a  film  type  of  cooling  is  obtained.  The  mass  flow 
of  water  into  the  chamber  amounts  to  less  than  1%  of  the  propellant  flow. 

The  total  pressure  probe  was  used  only  during  the  initial  steady  Btate  tests 
to  verify  the  design  Mach  number  at  the  nozzle  entrance. 

The  propellant  feed  system  is  pressure  fed  by  a  battery  of  trailer - 
mounted  dry  nitrogen  cylinders.  In  addition  to  tank  pressurization,  nitrogen 
gas  is  used  for  pneumatic -ope rated  valves,  post  firing  propellant  line  purg¬ 
ing  and  for  operation  of  pressure  receiver-transmitters.  The  nitrogen  is 
distributed  throughout  the  system  by  an  array  of  Grove  Power- reactor  dome 
type  pressure  regulators.  The  outlet  pressures  of  the  Grove  regulators 
are  controlled  by  the  dome  pressure  applied  to  them  through  individual 
venting-type  hand  loader  regulators  mounted  on  the  test  console. 


The  stainless- steel  liquid-oxygen  tank,  insulated  by  a  layer  of 
polystyrene,  is  filled  prior  to  test  from  a  low  pressure  vacuum -insulated 
250  gallon  storage  vessel.  The  fuel  tank  is  loaded  from  a  55  gallon  drum 
of  JP-5A.  In  addition  to  the  nitrogen  inlet  and  propellant  outlet,  each  tank 
is  equipped  with  electrically  operated  vent  valves,  pre-set  pressure  relief 
valves  and  tank  pressure  taps.  The  nitrogen  inlet  on  the  liquid  oxygen  tank 
is  fitted  with  a  diffuser  which  directs  the  pressurizing  medium  away  from 
the  surface  of  the  oxidizer.  This  is  done  to  prevent  agitation  and  subsequent 
boiloff  of  the  liquid  oxygen  and  to  minimize  the  likelihood  of.  the  diffusion 
and/or  dissolving  of  the  warm  nitrogen  gas  in  the  cryogenic  oxidizer. 

Propellant  flow  to  the  thrust  chamber  iB  controlled  through  individual 
propellant  valves  which  are  electrically  coupled  by  means  of  a  time  delay 
relay  circuit.  This  procedure  allows  for  control  of  the  lead  between  oxidizer 
and  fuel  injection.  Secondary  shut-off  valves  are  provided  in  each  of  the 
propellant  feed  lines  as  a  precaution  against  malfunction  of  the  main  pro¬ 
pellant  valves.  Cavitating  venturis  are  employed  in  the  fuel  and  oxidizer 
lines  to  regulate  the  propellant  flow  and  eliminate  feed  line  coupling  in  the 
event  of  chamber  pressure  oscillations.  A  special  igniter  has  been  designed 
(Fig.  3.  5)  to  insure  rapid  ignition  of  the  propellants  upon  injection  into  the 
combustion  chamber.  The  igniter  consists  of  a  paraffin  base  end  burning 
propellant  grain,  a  cupron  ignition  wire  and  a  copper  fuse  wire.  The  copper 
fuse  wire  which  passes  across  the  end  of  the  grain,  completes  a  relay  circuit 
which  holds  the  main  propellant  valves  closed  until  the  fuse  wire  is  burned 
through  by  the  igniter.  The  presence  of  a  pilot  flame  upon  injection  of 
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propellants  eliminates  the  possibility  of  forming  the  highly  explosive 
substance  which  results  from  cold  flow  mixing  of  liquid  oxygen  and  JP-5A. 

Immediately  prior  to  firing,  the  liquid  oxygen  lines  and  injector 
manifold  are  chilled  by  flushing  the  cryogenic  fluid  through  the  system 
until  liquid  flow  is  indicated.  At  this  point  the  firing  button  is  depressed 
which  transfers  test  control  to  a  programmed  sequential  timer.  The  fol¬ 
lowing  commands  are  then  automatically  initiated. 

1.  Open  secondary  shut-off  valves. 

2.  Initiate  ignition  which  in  turn  causes  main  propellant  valves 
to  open  in  predetermined  sequence. 

3.  Initiate  experiments  and  operate  remote  auxiliary  equipment 
such  as  cameras,  pumps,  etc. 

4.  Close  secondary  shut-off  valves. 

5.  Open  purge  valves  which  causes  gaseous  nitrogen  to  flush 
propellant  lines,  injector  and  chamber  of  residual  propellants. 

6.  Return  all  systems  to  pre-firing  condition. 

Instrumentation  requirements  are  divided  into  two  categories; 

monitoring  systems  and  permanent  recording  systems.  Monitoring  infor¬ 
mation,  such  as  tank  pressures,  injection  pressures,  chamber  pressure 
and  valve  position  is  transmitted  to  the  control  building  where  it  is  dis¬ 
played  on  the  test  console.  Those  data  that  require  permanent  recording 
such  as  propellant  flow  rates,  thrust,  chamber  pressure  gradients  and  the 
output  of  high  frequency  response  pressure  transducers  used  in  instability 
studies  are  transmitted  to  the  instrumentation  facility. 
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Monitoring  of  the  various  pressures  is  facilitated  with  the  use  of 


Ashcroft  pressure  receiver-transmitters.  These  systems  sense  the  desired 
pressure  (0-i000psig)  at  its  source  and  uses  the  information  to  control  the 
output  of  a  low  pressure  {0-15  psig)  regulator.  This  signal  in  turn  is  trans¬ 
mitted  to  the  test  console  where  it  is  displayed  on  a  companion  pre-calibrated 
Bourdon  gage  whose  face  markings  correspond  to  the  source  pressure  range. 
Thus,  several  hundred  feet  of  tubing  containing  either  liquid  oxygen,  JP-5A 
or  combustion  products  have  been  replaced  with  an  equivalent  tube  length  of 
low  pressure  gaseous  nitrogen. 

Instrumentation  for  acquiring  steady  state  axial  pressure  gradients, 
which  are  used  to  verify  theoretical  spray  combustion  models,  categorize 
various  injector  configurations  and  serve  as  boundary  conditions  for  wave 
propagation  studies,  was  selected  after  consideration  of  several  factors, 
namely;  proximity  of  measurements,  expected  pressure  gradient,  accuracy, 
and  cost.  Assuming  that  the  pressure  was  to  be  measured  at  eight  longi¬ 
tudinal  locations,  the  apparent  solution  would  be  to  install  eight  transducers 
down  the  length  of  the  chamber.  However,  each  of  the  considerations  above 
eliminated  this  technique. 

Testa  conducted  on  the  shorter  length  chambers  would  require  a 
pressure  measurement  at  one  inch  intervals.  The  physical  size  of  the 
diaphragm  on  a  water  cooled  transducer  varies  from  1/2  to  1  inch  depending 
on  the  manufacturer.  Therefore,  each  transducer  would  actually  measure  an 
average  pressure  over  an  interval  comparable  to  the  center  to  center  distance 
between  pressure  taps.  In  addition,  it  was  expected  that  the  gradient  at 
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certain  sections  of  the  motor  would  be  as  small  as  1  psi  in  200  psi 
between  adjacent  taps.  Using  a  500  psi  full  range  transducer,  a  drift 
of  .  1%  of  full  scale  in  opposite  directions  by  two  adjacent  transducers 
would  yield  zero  pressure  gradient.  Rather  than  develop  a  method  of 
compensating  for  the  different  drift  characteristics  of  adjacent  transducers 
and  associated  electronics,  it  appeared  more  feasible  to  develop  a  tech¬ 
nique  that  employed  a  single  transducer  to  measure  the  pressure  gradient. 

In  addition,  if  such  a  system  could  be  developed  the  dollar  savings  to  the 
experimental  program  would  be  significant.  Rather  than  employ  eight  water 
cooled  transducers  it  would  be  possible  to  use  a  single  transducer.  It  will 
subsequently  be  shown  that  the  transducer  need  not  be  cooled,  resulting  in 
an  additional  savings.  All  of  the  above-mentioned  factors  lead  to  the  devel¬ 
opment  of  a  pressure  scanner  (Fig.  3.  6). 

The  scanner  permits  a  single  transducer  to  travel  from  pressure 
tap  to  pressure  tap  and  consecutively  record  the  pressure  at  eight  different 
locations  in  the  chamber.  The  pressure  scanner  consists  of  three  piston 
and  cylinder  units.  Units  A  and  B  contain  facilities  for  measuring  the  pres¬ 
sure  at  eight  locations,  while  ui.it  C  is  a  hydraulic  driver.  The  taps  located 
in  the  rocket  motor  are  extended  to  the  bulkhead  of  either  unit  A  or  B  and 
then  to  compartments  formed  by  the  inner  cylinder  or  piston,  "o"  rings,  and 
the  outer  travelling  cylinder  which  houses  a  transducer.  The  reciprocating 
motion  of  the  outer  cylinders,  (which  enables  the  transducer  to  consecutively 
measure  the  pressure  at  adjacent  taps  in  the  rocket  motor)  is  governed  by  a 
double  acting  hydraulic  piston  and  cylinder;  unit  C.  Provisions  have  been 
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made  for  varying  the  speed  and  length  of  stroke  of  the  unit,  by  varying  the 
hydraulic  pressure  and  position  of  the  reversing  rnic roswitche s  which  in 
turn  controls  the  hydraulic  fluid  feed  and  vent  valves.  In  addition,  a  method 
of  continuously  recording  the  exact  location  of  the  transducer  housing  has 
been  developed.  This  consists  of  measuring  the  voltage  drop  between  a 
moving  electrical  contact  mounted  on  the  transducer  housing  and  a  flat 
copper  bar  mounted  alongside  the  scanner  cylinder.  The  copper  bar  is 
actually  an  eight  stepped  voltage  divider.  The  length  of  each  step  is  equal 
to  the  length  of  the  compartments  in  the  scanner  cylinder.  The  distance 
between  steps  or  sections  is  equal  to  the  thickness  of  the  "O"  rings  sepa¬ 
rating  the  compartments.  Each  of  the  copper  bar  sections  is  at  a  different 
predetermined  voltage.  Therefore,  by  simultaneously  recording  the  output 
of  the  transducer  and  the  voltage  drop  to  the  moving  electrical  contact,  it 
is  possible  to  match  the  pressure  magnitude  with  the  location  at  which  it 
was  measured. 

The  scanner  was  subsequently  modified  by  replacing  the  hydraulic 
driving  force  with  a  variable  speed  AC  motor  coupled  to  a  slider  crank 
mechanism,  (Fig.  3.  7).  More  recently  a  pressure  scanning  system,  manu¬ 
factured  by  the  Scanivalve  Co.  ,  San  Diego,  California  has  been  installed. 

The  unit  is  basically  a  rotary  version  of  the  previous  scanners,  which 
couples  a  single  transducer  to  as  many  as  48  pressure  sources  through  a 
rotating  port. 

The  output  of  the  transducers  is  recorded  on  a  recording  oscillo¬ 
graph.  To  record  the  absolute  pressure  at  each  port  would  result  in 
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reduced  resolution.  Therefore,  only  the  pressure  difference  about 
a  fixed  point  is  recorded.  Thus  40  psi  rather  than  200  psi  is  Till-scale 
deflection  on  the  recorder.  The  above  technique  is  accomplished  by  sev¬ 
eral  methods  depending  upon  the  type  of  transducer  employed  in  the  scan¬ 
ner.  Using  a  piezoelectric,  crystal  transducer,  the  crystal  is  grounded 
until  the  starting  transient  nas  been  completed.  The  pressure  on  the 
transducer  diaphragm  at  the  instant  that  the  ground  connc  .  tion  is  broken 
corresponds  to  zero  output.  Therefore,  the  subsequent  output  of  the 
transducer  is  proportional  to  the  difference  between  the  pressure  being 
measured  and  the  pressure  on  the  transducer  when  the  ground  connection 
was  broken.  To  determine  the  absolute  pressure  at  each  port,  from  the 
records  of  the  gradients,  requires  a  measurement  of  the  absolute  pressure 
at  one  point  in  the  chamber. 

Experimental  data  on  wave  shape  behavior  is  obtained  with  three 
high  frequency- response,  flush-mounted,  water-cooled  pressure  trans¬ 
ducers  (Photocon  model  no.  352)  spaced  axially  along  the  chamber.  Data 
are  recorded  on  a  frequency  modulated  magnetic  tape  recorder- reproducer 
at  60  ips.  The  data  is  played  at  1  ips  into  a  recording  oscillograph  with 
a  paper  speed  of  25  ips.  The  final  data  record  thus  has  an  expanded  time 
scale  of  1500  ips.  Frequency- spectrum  analysis  of  the  waves  can  be 
accomplished  by  connecting  the  output  of  the  tape  recorder  playback 
amplifiers  to  a  Panoramic  Model  LP-la  audio-frequency  spectrum  analyzer. 

A  continuous  record  is  made  of  fuel  and  oxidizer  flow  rates  during 
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the  test  firing.  A  turbine-type  flowmeter  is  installed  in  each  of  the 
propellant  feed  lines.  The  output  of  the  flow  transducers  is  an  AC 
signal  whose  frequency  is  proportional  to  the  propellant  flow  rate. 

The  sinusoidal  signal  is  converted  to  a  pulse  signal  of  equal  frequency 
which  is  then  applied  through  a  stylus  to  an  electrosensitive  paper 
recorder.  A  60  cps  calibration  signal  is  simutaneously  recorded  for 
reference  purposes. 


SECTION  IV 


EXPERIMENTAL  RESULTS 

A.  Scope  of  the  Program 

An  experimental  program  v'<.  s  conducted  in  order  to  obtain  steady 
state  and  v/ave  propagation  data  in  a  liquid  propellant  rocket  syeten  Steady 
state  experiments  serve  to  examine  the  validity  of  spray  combustion  models 
as  w  1  as  provide  a  source  for  empirical  data  on  droplet  size  distributions. 

In  addition,  boundary  and  initial  conditions  for  theoretical  and  experimental 
wave  propagation  studies  are  obtained.  The  wave  propagation  experiments 
are  for  the  purpose  of  determing  the  parameters  which  influence  energy 
oupling  to  the  wave,  wave  coalescence,  and  the  effect  of  velocity  and  pres¬ 
sure  gradients  on  wave  attenuation. 

B.  Steady  State 

In  the  design  of  liquid  propellant  rocket  combustors  it  is  desired 

to  have  complete  burning  take  place  within  the  combustion  chamber.  In 

addition,  the  liquid  and  gasdynamic  flow  fields  should  be  inherently  stable. 

4  34 

Recent  investigations  ’  and  Section  II  of  this  report,  resulted  in  methods 
to  determine  the  minimum  chamber  length  required  to  insure  complete 
combustion  and  the  flow  fields  response  to  input  disturbances.  The  theo¬ 
retical  analysis  predicted  mass  liberation  schedules,  gs^s  and  droplet  velocity 
histories,  and  axial  pressure  histories.  It  was  found  that  increases  in  the 
mass  liberation  rate  are  accompanied  by  attendant  increases  in  the  magnitude 
of  the  pressure  gradient  with  a  corresponding  increase  in  the  velocity  gradient. 
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Initial  experiments  designed  to  measure  axial  pressure  gradients  were 

4 

in  good  agreement  with  the  theory  for  a  given  set  of  injection  parameters,  . 

In  order  to  obtain  correlation,  it  was  necessary  to  assume  a  drop  size 
distribution  and  a  droplet  formation  distance,  i.  e.  ,  the  required  distance 
from  the  injector  face  for  the  injected  ligaments  to  form  a  discrete  droplet 
distribution.  The  analysis  distinctly  indicated  that  energy  (mass)  liberation 
gradients  can  be  deduced  from  pressure  gradient  measurements.  Thus, 
the  chamber  length  for  optimum  propellant  utilization  and  thrust  production 
is  related  to  the  axial  pressure  gradient.  In  a  motor  whose  geometry  and 
propellant  mass  flow  rate  are  fixed,  this  gradient  is  almost  entirely  depen¬ 
dent  on  the  nature  of  the  injection  process. 

In  addition,  the  behavior  of  an  input  pressure  disturbance  is  depen¬ 
dent  on  the  gasdynamic  velocity  field  through  wh’ ch  the  disturbance  propagates. 
Thus,  the  characteristic  gradients  or  signatures  produced  by  various  in¬ 
jectors  are  important  parameters  in  determining  the  stability  or  inst¬ 
ability  of  a  thrust  chamber  system.  Therefore,  an  experimental  program 
was  initiated  in  order  to  determine  the  pressure  gradient  characteristics 
of  various  injector  configurations. 

The  500  lbf  nominal  thrust  rocket  engine  system  described  in 
Section  III  wae  used  for  this  investigation.  The  chamber  length  was  varied 
between  8  and  24  in.  ,  although  all  of  the  data  reported  here  are  for  17  in. 
chambers  except  where  noted.  The  chamber  and  nozzle  throat  diameters 
are  2  and  1. 65  in.  ,  respectively,  resulting  in  a  contraction  ratio  of  1. 47. 

The  propellants  are  JP-5A  and  liquid  oxygen. 
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All  of  the  injectors  used  in  this  investigation  contain  16  fuel  and 
16  oxidizer  orifices.  The  oxidizer  orifice  is  maintained  constant  at 
0.  0595  in.  diameter,  and  the  fuel  orifices  are  0.  042  in.  diameter,  except 
in  the  S3  showerhead  injector,  which  has  m  .ltisized  orifices  in  order  to 
obtain  a  greater  distribution  of  fuel  droplet  sizes.  All  impinging  injectors 
are  of  the  unlike  impinging  stream  pattern  type,  i.  e.  ,  fuel  on  oxidizer. 
Drawings  of  the  various  injector  configuration  are  presented  in  Fig.  4.  1. 

A  summary  of  the  above  injectors  is  tabulated  in  Fig.  4.  2. 

In  order  to  eliminate  the  effects  of  injection  velocity  on  the  drc 
ballistics  and  subsequent  evaporation  history,  the  total  liquid  propc 
mass  flow,  oxidizer-fuel  ratio,  and  total  injection  area  of  both  Jl 
liquid  oxygen  are  held  constant  throughout  the  series  of  tests  present 
reported.  JP-5A.  injection  velocity  is  85  fps.  The  average  total  propellant 
flow  is  2.  35  lbm/sec  at  an  0/F  ratio  of  2.  62.  A  5%  variation  from  the 
average  flow  rate  is  allowed  for  inclusion  in  the  test  results. 

The  data  for  the  injectors  tested  are  presented  both  as  pressure 
vs.  axial  distance  and  pressure  gradient  (AP/6x)  vs.  the  average  axial 
distance  over  which  the  gradient  was  computed.  It  should  be  noted  that 
later  testa  included  a  static  pressure  tap  at  the  injector  face.  Consequently, 
the  more  recent  data  presented  includes  a  pressure  measurement  at  this 
point. 

Figure  4.  3  shows  the  data  for  two  showerhead  injectors  and  one 
impinging  injector.  The  two  showerheads  both  contain  four  alternating 
rings  of  fuel  and  oxidizer,  but  differ  in  fuel  orifice  diameter.  The  S2 
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injector  contains  sixteen  0.  042-in.  -diameter  orifices  whereas  the  S3 
has  eight  0.  035-in.  -,  four  0.  043-in.  -,  and  four  0.  052-in.  -diameter 
orifices.  The  impinging  injector  was  designed  for  an  impingement 
distance  (x^ )  of  0.  28  in.  The  data  indicate  a  much  steeper  gradient 
and,  hence,  a  more  localized  region  of  evaporation  and  combustion  for 
the  impinging  injector  than  for  either  of  the  showerheads.  The  use  of 
multisized  orifices  results  in  a  more  uniform  pressure  gradient  and  hence 
a  spreading  of  the  evaporation  and  combustion  zone. 

Data  for  a  radial  sheet  and  tangential  sheet  impinging  injector  with 
the  same  design  impingement  distance  are  shown  in  Fig.  4.4.  Both  injectors 
contain  the  same  number  and  3ize  of  orifices.  The  tangential  sheet  injec¬ 
tor  contains  a  single  ring  of  doublets,  whereas  the  radial  sheet  injector 
contains  two  rings  of  doublets.  The  latter  injector  increases  mixing  and 
atomization  of  the  fuel  and  oxidizer  and  results  in  more  rapid  utilization 
of  the  propellants.  Thus,  the  maximum  point  on  the  pressure  gradient 
curve  for  the  radial  sheet  injector  is  closer  to  the  injector  face  than  the 
corresponding  point  for  the  tangential  sheet  injector.  This  indicates  that 
the  main  evaporation  and  combustion  zone  has  been  shifted  upstream  with 
the  12  injector  and  implies  that  a  shorter  chamber  would  be  required  to 
obtain  maximum  propellant  evaporation. 

The  results  obtained  with  two  additional  impinging  injectors,  each 
of  which  contains  three  impingement  distances,  are  shown  in  Fig.  4.  5  to¬ 
gether  with  the  single  impingement  point,  16  injector.  As  the  impingement 
point  is  moved  downstream,  the  degree  of  mixing  decreases  and  the 


evaporation  and  combustion  zone  is  spread  out  over  a  greater  axial 
distance.  The  high  amplitude  and  thin  width  of  the  pressure  gradient 
peak  for  the  16  injector  indicate  a  zone  of  high  evaporation  and  combus¬ 
tion  approximately  4  in.  from  the  injector  face.  The  injector  with 
impingement  distances  less  than  1  in.  has  a  peak  in  the  pressure  gradient 
curve  6  in.  from  the  injector  face,  whereas  the  pressure  gradient  for  the 
last  injector  indicates  a  peak  near  the  nozzle  entrance.  This  implies  a 
high  degree  of  evaporation  and  gas  generation  in  the  downstrean  portions 
of  the  chamber. 

Data  have  been  obtained  from  three  tangential  sheet  impinging 
injectors  with  design  impingement  points  greater  than  1  in.  ,  Fig.  4.  6. 
Although  it  cannot  be  readily  determined  whether  or  not  impingement 
takes  place  at  such  large  distances  from  the  face,  the  data  indicate  that, 
as  the  intended  impingement  point  is  moved  downstream,  the  rate  of  pres¬ 
sure  drop  decreases.  At  the  maximum  impingement  distance  tested,  the 
effect  of  poorer  mixing  and  low  propellant  utilization  is  evidenced  by  the 
decrease  in  over-all  chamber  pressure,  despite  a  somewhat  higher  propel¬ 
lant  flow  rate. 

The  performance  with  this  last  injector  (19)  is  inferior  to  that 
obtained  with  the  showerhead  S2  injector  discussed  previously.  Aside  from 
the  obvious  differences  between  the  two  injectors,  the  showerhead  contains 
four  alternating  rings  of  fuel  and  oxidizer,  whereas  the  impinging  injector 
contains  a  single  ring  of  doublets.  It  appears  then  that  distribution  of  the 
propellant  over  a  greater  portion  of  the  chamber  cross  section  enhances 
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mixing  and  atomization  and  promotes  the  evaporation  process.  Test6 
conducted  with  a  two-ring  showerhead  indicate  that  a  chamber  length  of 
21  in.  is  required  for  optimum  performance,  whereas  the  four-ring  injec¬ 
tor  only  required  14  in. 

The  three  injectors  just  discussed  are  equipped  with  pressure 
taps  in  the  injector  face.  The  data  show  that  the  point  of  maximum  chamber 
pressure  is  downstream  from  the  injector  face,  and  the  pressure  gradient 
near  the  face  is  positive.  This  is  indicated  by  a  negative  value  of  -(AP/Ax). 
This  suggests  the  existence  of  a  recirculation  zone  or  negative  velocity 
gradient  in  the  region  near  the  injector  face. 

Variations  in  injector  configuration  produce  variations  in  axial  pres¬ 
sure  history  in  a  constant  area  liquid  propellant  rochet  motor.  These 
pressure  variations  are  due  to  changes  in  the  mixing  and  shattering  char¬ 
acteristics  of  different  injectors  and  are  concurrent  with  unique  evaporation 
and  mass  liberation  profiles.  For  an  injector  that  promotes  rapid  droplet 
breakup  and  intimate  mixing  of  propellants,  such  as  the  16  used  in  these 
experiments,  the  chamber  pressure  has  a  maximum  gradient  close  to  the 
injector  face  with  a  correspondingly  high  gas  acceleration.  This  action 
exerts  additional  gasdynamic  effects  on  the  unconsumed  droplets,  enhancing 
propellant  consumption. 

Since  the  attenuation  or  amplification  of  a  pressure  wave  is  depen¬ 
dent  on  the  gasdynamic  field  through  which  it  propagates,  the  inherent  sta¬ 
bility  of  an  injector -chamber  system  is  keyed  to  the  pres  sure -velocity 
gradient  produced. 
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In  addition,  the  data  can  be  used  to  establish  design  criteria  for 
steady-state  operation  and  to  predict  critical  conditions  that  can  exist 
because  of  large  pressure  gradients  that  result  in  high  heat-transfer  rates 
to  the  injector  face  or  rocket  chamber  wall. 

C.  Wave  Propagation  Phenomena 

Longitudinal  wave  propagation  studies  were  conducted  in  a  liquid 
propellant  rocket  motor  in  order  to  observe  wave  deformation  phenomena 
and  to  define  the  parameters  that  determine  whether  an  input  disturbance 
will  attenuate  or  amplify.  Since  combustion  instability  is  a  measure  of 
the  relative  amounts  of  energy  accumulation  in  a  cavity  in  contradistinction 
to  the  energy  dissipation  from  the  cavity,  the  mechanisms  that  allow  such 
behavior  should  be  analyzed  in  detail.  Particular  emphasis  must  be  placed 
upon  the  interaction  of  pressure  waves  and  the  fluid  dynamic  field  in  ac¬ 
cordance  with  remarks  made  in  the  introduction,  which  are  repeated  here 
for  completeness. 

Energy  coupling  between  the  propagating  wave  and  the  combustion 
process  is  dependent  upon  the  relative  magnitudes  of  characteristic  times 
associated  with  the  source  and  sink.  If  the  relaxation  time  of  a  significant 
transport  process  in  the  conversion  of  liquid  propellants  to  gaseous  products 
is  less  than  or  equal  to  the  wave  residence  time  in  a  volume  element, 
coupling  can  occur  with  resultant  unstable  operation.  The  leading  edge 
of  a  passing  wave  in  a  reacting  droplet  system  can  increase  the  rate  of 
evaporation,  which  couples  as  a  mass  source  to  the  trailing  edge  of  the 
passing  wave,  thus  producing  amplification.  It  can  further  be  seen  that, 


for  a  long  relaxation  time,  the  mass  source  can  generate  wavelets  that 
coalesce  as  they  propagate  and  ultimately  overtake  the  initial  wave  that 
caused  the  disturbance,  again  causing  amplification.  In  addition  as  a  wave 
propagates  in  a  gas,  it  deforms.  The  manner  and  extent  of  the  deformation 
is  dependent  on  the  type  of  wave,  compression  or  expansion,  and  the  velocity 
gradient  of  the  gas  through  which  it  propagates.  Therefore,  as  a  wave  moves 
in  a  rocket  motor  it  changes  its  geometry  and  residence  time  in  a  volume 
element,  which  in  turn  affect  the  nature  of  the  energy  or  mass  coupling  to 
the  propagating  wave  and  the  ultimate  stability  of  the  system. 

Initial  wave  propagation  experiments  were  conducted  in  a  2-in.  - 
diameter,  500-lbf  nominal  thrust,  JP-5A  liquid  oxygen  rocket  motor.  The 
injector  is  of  the  shower-head  type  with  16  fuel  and  16  oxidizer  orifices, 
similar  to  the  S2  injector  described  previously.  The  data  discussed  below 
is  obtained  with  <.  total  propellant  flow  of  2.  37  lbm/ sec  at  an  oxidizer-fuel 
ratio  of  2.  79.  The  chamber  length  is  22.  5  in.  measured  from  the  injector 
face  to  the  start  of  the  converging  nozzle.  The  nozzle  contraction  ratio  is 
1.47;  resulting  in  a  high  Mach  number  profile  through  the  chamber.  High 
frequency- response  pressure  transducers  are  flush-mounted  in  the  chamber 
at  3,  13,  and  21  in.  from  the  injector  face.  Approximately  six  seconds 
after  startup  (steady-state  chamber  pressure  is  165  psia  measured  3  in. 
from  the  injector  face),  the  wave  generator  tube  diaphragm  is  ruptured, 
and  a  wave  propagates  into  the  rocket  motor  thrust  chamber.  The  wave 
generator  is  actually  a  modified  shock  driver  tube  mounted  axially  through 
the  injector  manifold  and  utilizing  the  rocket  motor  combustor  as  the  driven 


tube.  The  transition  piece  between  the  shock  driver  tube  and  the  combustor 
shapes  the  input  wave.  The  diaphragm  is  in  a  plane  parallel  to  the  injector 
face.  Diaphragm  rupture  is  accomplished  with  a  solenoid  actuated  needle. 

The  wave  generator  tube  driver  pressure  is  725  psia. 

Initially,  the  diaphragm  was  mounted  within  one  quarter  of  an  inch 
of  the  injector  face  (Fig.  4.7).  The  diaphragms  were  fabricated  from  stain¬ 
less  steel  or  brass  shim  stock  in  a  variety  of  thickness  ranging  from  0.  005  in. 
to  0.  020  in.  The  close  proximity  to  the  combustion  gases  often  caused 
failure  of  the  diaphragm  during  the  rocket  motor  starting  transient.  At 
best,  the  high  thermal  load  imposed  resulted  in  a  severely  stressed  dia¬ 
phragm  which  did  not  petal  upon  bursting.  As  a  result  the  input  wave  was 
not  reproducible.  Modifications  to  the  system  resulted  in  moving  the  dia¬ 
phragm  upstream  of  the  injector,  and  using  a  transition  piece  between  the 
driver  tube  and  injector  face  to  shape  the  wave.  A  schematic  of  the  thrust 
chamber-wave  generator  system  is  shown  in  Fig.  4.8.  A  typical  injector, 
modified  to  accept  the  wave  generator  transition  piece  is  shown  in  Fig.  4.  9. 

The  pressure  time  data  of  a  typical  run  is  shown  in  Fig.  4.  10.  The 
three  traces  correspond  to  pressure  recordings  taken  3,  13,  and  21  inches 
from  the  injector  face.  The  data  separates  into  two  parts.  One  part  con¬ 
sists  of  the  input  wave  from  the  wave  generator  which  oscillates  axially 
in  the  chamber.  This  is  shown  in  the  upper  part  of  the  figure.  The  second 
part  is  a  spontaneous  oscillation  which  occurs  soon  after  the  first  wave  group 
attenuates.  This  is  shown  in  the  lower  part  of  the  picture.  These  latter 
waves  are  discussed  first. 
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The  output  of  all  transducers  indicates  a  periodic  disturbance 
propagating  longitudinally  in  the  thrust  chamber  (Fig.  4.  10b).  The  wave 
is  initiated  at  the  injector  face,  propagates  downstream,  reflects  from  the 
sonic  plane  of  the  nozzle  throat,  and  propagates  upstream.  At  the  injector 
face,  the  wave  is  reflected  once  more  and  travels  toward  the  nozzle.  The 
period  of  the  disturbance  (defined  as  the  time  between  two  successive  waves 
traveling  in  the  same  direction)  is  1.  37  msec.  The  amplitude  of  the  over¬ 
pressure  is  approximately  100  psi  with  a  rise  time  of  0.  05  msec.  The  pres¬ 
sure  decays  behind  each  incident  wave  to  the  steady- state  chamber  pressure 
before  the  time  of  arrival  of  the  reflected  wave.  The  frequency  of  the  dis¬ 
turbance  is  730  cpsr  With  a  chamber  length,  measured  from  the  injector 
face  to  the  sonic  plane,  equal  to  2  ft,  the  equivalent  wavelength  would  appear 
to  be  4  ft.  Therdore,  the  expected  wave  velocity  relative  to  the  gas  would  be 
2920  fps.  This  is  below  the  expected  speed  of  sound  for  liquid  oxygen-JP-5A 
combustion  products  and  is  less  than  the  propagation  velocity  as  determined 
from  the  transducer  outputs  in  the  following  manner. 

Velocities  of  wave  propagation  relative  to  a  stationary  observer  are 
determined  by  dividing  the  distance  between  adjacent  transducers  by  the  time 
required  for  the  wave  to  travel  between  transducers.  The  average  value  of 
the  upstream  and  downstream  velocities  of  propagation  yields  the  propaga¬ 
tion  velocity  relative  to  the  gas  which,  in  this  case,  is  3350  fps.  Previous 
4 

analytical  work  has  shown  that  the  droplet  and  gas  ballistics  are  consistent 
with  a  chamber  temperature  of  5500° R  and  a  sonic  velocity  of  3580  fps.  Since 
the  temperature  and  speed  of  sound  are  decreasing  through  the  converging 


portion  of  the  nozzle,  the  average  sonic  velocity  relative  to  the  gas 
would  be  less  than  3580  fps. 

The  apparent  anomaly  between  the  wave  velocity  as  determined 

from  the  frequency  and  equivalent  wavelength  (2920  fps)  and  the  velocity 

as  determined  from  the  upstream  and  downstream  averaging  (3350  fps) 

can  be  explained  by  work  done  pieviously  at  the  Polytechnic  Institute  of 
35 

Brooldyn,  .  It  was  shown  that  the  resonant  frequency  for  a  duct  contain¬ 
ing  a  gas  flow  at  some  finite  Mach  number  M  is  given  by 

f  =c(l  -  M2)/2£ 
r 

or  the  equivalent  wavelength  is 

X  =  [2A/(1  -  M3)] 

where  j i  is  the  geometrical  length  of  the  duct,  c  is  the  average  wave 
velocity  realtive  to  the  gas,  and  the  duct  is  assumed  to  behave  as  a  half¬ 
wave  resonator.  Therefore,  the  product  of  the  measured  frequency  and 
twice  the  duct  length  does  not  yield  the  average  wave  velocity,  but  rather 
this  latter  quantity  multiplied  by  a  factor  of  (1  -  M3).  For  a  resonant  fre¬ 
quency  of  730  cps,  a  duct  length  of  2  ft,  and  an  average  wave  velocity  of 
3350  fps,  the  Mach  number  as  determined  from  the  previous  equations  is 
0.  367.  Since  the  nozzle  entrance  Mach  number  on  the  motors  used  in  these 
tests  is  0.  45,  and  the  duct  length  includes  the  convergent  portion  of  the 
nozzle,  the  average  Mach  number  as  determined  from  the  resonant  frequency 
appears  to  be  reasonable. 

The  method  described  above  for  calculating  average  wave  velocities 
is  also  used  to  calculate  average  gas  velocities  which  in  turn  yields  the 
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average  Mach  number,  These  calculations  are  shown  in  Fig.  4.  11.  This 
latter  Mach  number  was  compared  to  values  obtained  above  from  frequency 
measurements  and  it  is  seen  that  the  values  agree  very  well,  (e.  g.  ,  0.  367 
vs.  0.  363).  .However,  these  localized  measurements  bring  out  another 
significant  feature.  If  the  local  sonic  velocity  is  used  to  indicate  the  local 
gas  temperature,  then,  during  this  type  of  oscillation,  it  is  seen  that  the 
gases  in  the  injector  end  of  the  motor  are  cooler  than  those  at  the  nozzle 
end.  In  other  words,  there  is  a  variable  stagnation  temperature  in  the 
chamber  due  to  different  evaporation  rates  of  the  propellants  and  consequent 
variable  0/F  ratio  during  combustion.  In  any  detailed  wave  analysis,  both 
the  fuel  and  droplet  ballistics  must  be  considered  to  give  accurate  wave 
behavior. 

It  is  obvious  that  for  low  Mach  numbers  and  large  contraction  ratios 
the  effect  on  the  resonant  frequency  is  small,  but  as  the  Mach  number 
increases,  the  reduction  in  natural  or  resonant  frequency  increases.  If 
one  has  a  beforehand  knowledge  of  the  wave  velocity  relative  to  the  gas, 
say  3350  fps,  and  uses  simple  acoustic  theory  (M  =  0),  the  resonant  fre¬ 
quency  for  a  2-ft  chamber  would  be  840  cps.  As  the  Mach  number  is  increased 
to  0. 367,  a  13%  decrease  in  resonant  frequency  is  noted. 

An  important  conclusion  to  be  derived  from  the  foregoing  discussion 
is  that  acoustic  theories  will  not  give  accurate  quantitative  results  of  the 
analysis  of  w ave  oscillations  in  rocket  motors,  especially  those  having  small 
contraction  ratios  with  resultant  high  Mach  number  profiles. 

Now  consider  the  initial  part  of  the  wave  train,  (Fig.  4.  10a).  A 

65 


J 


diagram  showing  the  time  of  arrival  of  the  wave  at  the  various  transducer 
locations  is  shown  in  Fig.  4.  12,  along  with  the  tabulation  of  various  para¬ 
meters.  The  initial  wave  is  generated  with  a  wave  slope  of  32  psi/in.  A6 
the  wave  propagates  through  the  combustion  zone,  i.  e.  ,  between  3  and  10  in. 
from  the  injector  face,  where  large  longitudinal  pressure  and  velocity 
gradients  are  present,  it  is  seen  to  broaden.  Actual  transducer  outputu 
from  the  3-  and  13-in.  locations  reveal  a  train  of  secondary  wavelets  fol¬ 
lowing  behind  the  initial  input  wave  (Fig.  4.  10a).  The  pressure  history 
at  the  21 -in.  transducer  indicates  that  all  of  the  secondary  waves  have 
coalesced  with  the  initial  wave  to  forma  single  steepening  wave.  The  sec¬ 
ondary  wavelets  have  a  frequency  below  that  expected  for  a  radial  mode  and 
are  not  characteristic  of  tangential  modes.  They  are  presumed  to  be  due  to 
the  interaction  of  the  axial  wave  with  the  evaporation  zone,  resulting  in  a 
shattering  of  drops  with  an  attendant  increase  in  over-all  evaporation  and 
combustion  rates.  Since  the  steady- state  worV  has  previously  shown  the 
evaporation  and  combustion  zone  to  be  limited  to  the  first  10  to  12  in.  of 
chamber  for  the  injector  currently  in  use,  it  would  not  be  expected  that 
combustion  waves  be  produced  in  the  latter  half  of  the  chamber.  It  is  inter¬ 
esting  to  note,  however,  that  the  secondary  or  combustion  waves  produced 
in  the  upstream  portion  of  the  chamber  coalesce  into  the  incident  wave  before 
arrival  at  the  21 -in.  transducer. 

The  wave  that  is  reflected  upstream  from  the  nozzle  end  is  initially 
attenuated.  However,  as  it  propagates  upstream,  it  steepens  considerably 
(from  17  to  105  psi/in.  ).  An  extrapolation  of  the  wave  diagram  for  the 
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incident  and  reflected  wave  indicates  that  the  equivalent  reflecting  sur¬ 
face  is  exactly  coincident  with  the  sonic  plane  at  the  nozzle  throat.  Further 
analysis  of  the  reflected  wave  reveals  that  the  overpressure,  (AP),  increases 
from  60  psi  to  193  psi.  This  is  measured  from  a  base  pressure  of  about 
310  psia.  Thus,  the  new  "steady  state"  pressure  increased  fron.  165  to 
310  psia  and  the  overpressure  from  60  to  193  or  by  130  psi.  Put  another 
way,  the  instantaneous  chamber  pressure  increased  from  165  to  503  psia. 
This  behavior  can  only  come  about  by  increased  mass/energy  accumulation 
in  the  rocket  combustion  chamber.  Thus,  the  experiments  indicate  that 
increased  propellant  evaporation  and  subsequent  burning  is  the  mechanism 
by  which  the  gasdynamic  flow -field  drives  the  waves. 

At  the  injector  face,  the  wave  is  reflected  and  starts  a  second 
traverse  of  the  chamber.  In  all,  three  complete  cycles  are  required  before 
the  steady-state  component  of  the  pressure  is  reduced  to  the  prewave  value. 

'  nother  item  of  interest  is  the  following.  The  local  speed  of  sound 
of  the  initial  wave  referenced  to  the  gas  is  given  in  column  1  of  Fig.  4. 13. 
The  data  presented  was  obtained  from  four  separate  tests,  the  first  three 
of  which  resulted  in  spontaneous  instability  (Fig.  4.  10b).  The  wave  is 
propagating  into  a  highly  oxidant  rich  combustion  gas  whose  temperature  is 
significantly  below  the  adiabatic  flame  temperature  for  the  injected  0/F 
ratio.  The  wave  velocity  progressively  increases  as  it  proceeds  downstream 
and  this  behavior  is  expected  as  a  result  of  increased  temperature  in  the 
axial  direction.  However,  and  more  important,  the  wave  continues  to  accel¬ 
erate  as  it  reflects  from  the  sonic  plane  of  the  nozzle  and  propagates  toward 
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the  injector.  This  indicates  that  the  temperature  of  the  gases  in  the  region 
near  the  injector  face  has  increased.  It  is  felt  that  this  behavior  can  be 
explained  by  the  increased  turbulence  of  the  gases  behind  the  wave  which 
induced  increased  mixing,  and  subsequent  combustion.  This  is  based  on 
the  fact  that  the  temperature  increase  produced  by  the  overpressure  due 
to  the  compressional  effect  of  the  wave  is  about  3%  and  this  value  is  too 
small  to  account  for  the  increased  wave  velocity.  In  test  4,  no  instability 
occurred  and  no  increase  in  velocity  was  observed. 

Now  consider  the  existence  of  possible  relaxation  times.  It  is  noted 
in  Fig.  4.  10  that  a  characteristic  time  of  about  80  microseconds  appears 
both  in  the  wave  produced  by  the  wave  generator  and  in  the  spontaneously 
generated  wave.  This  time  corresponds  to  the  period  of  the  first  tangential 
mode  of  oscillation.  If,  however,  the  first  tangential  mode  was  excited, 
then  it  would  propagate  axially  with  the  particle  velocity.  This  is  clearly 
not  the  case  in  both  wave  groups.  These  waves  propagate  at  the  local  speed 
of  sound.  Cold  tests  were  made  on  the  chamber  with  instrumentation  for 
longitudinal  and  tangential  and  radial  oscillatory  behavior.  These  data  did 
not  contain  this  significant  time  value.  The  pressure  transducer  and  their 
mounts  were  investigated  to  determine  if  the  oscillations  could  be  attributed 
to  a  mechanical  or  electrical  origin.  The  transducers  were  dynamically 
calibrated  in  a  shock  tube  both  for  amplitude  response  and  the  existence  of 
"ring".  In  addition,  hot  firings  were  conducted  with  transducers  mounted 
in  blind  holes.  The  results  indicated  that  the  80  microsecond  signal  was 
not  due  to  electrical  noise,  a  hard  mount,  or  a  faulty  transducer,  but  rather 
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occurs  in  the  combustion  gases.  Since  this  time  is  greater  than  the 
expected  chemical  relaxation  times  and  less  than  the  normal  sized  drop¬ 
let  diffusion  times  through  the  film  boundaries,  the  origin  of  the  signal 
is  in  question.  It  is  suggested  however,  that  a  non  steady  droplet  evapora¬ 
tion-diffusion  analysis  be  initiated  to  determine  fuel  vapor  diffusion  times 
of  various  sized  droplets  in  films  of  various  concentrations. 

The  previous  wave  propagation  experiments  were  conducted  pri¬ 
marily  for  observation  of  the  phenomena  and  to  deduce  any  mechanisms 
that  could  couple  to  the  wave  and  cause  sustained  oscillations.  The  next 
phase  of  the  program  was  to  conduct  a  systematic  investigation  to  determine 
the  conditions  under  which  input  disturbances  attenuate  or  amplify  in  liquid 
propellant  rocket  motors. 

Previous  work  has  demonstrated  the  importance  of  the  physical 
processes  of  atomization,  mixing  and  evaporation,  droplet  heat  transfer 
and  chemical  reaction  in  determining  the  stability  of  a  thrust  chamber.  If 
any  of  the  previous  processes  are  sensitive  to  pressure  waves,  then  it  is  pos¬ 
sible  that  coupling  occurs  with  the  result  that  oscillations  are  sustained  and 
amplified.  These  processes  are  not  distributed  throughout  the  entire  chamber, 
but  are  localized  in  extent.  However,  the  longitudinal  pressure  waves  pro¬ 
pagate  between  the  injector  face  and  the  sonic  plane  of  the  nozzle.  Since  it 
is  possible  for  the  wave  to  be  distorted  in  regions  of  the  chamber  containing 
little  or  no  sensitive  processes,  a  complete  understanding  of  wave  attenua¬ 
tion  phenomena  requires  particular  emphasis  upon  the  interaction  of  pressure 
waves  and  the  entire  thrust  motor  fluid  dynamic  field. 


The  experiments  reported  in  this  section  were  conducted  in  2  in. 
diameter  500  Ibf.  nominal  thrust,  JP-5A-LOX  rocket  motor.  The  injectors 
are  of  the  showerhead  type  or  unlike  impinging  type  with  16  fuel,  (0.  042-in. 
diam)  and  16  oxidizer,  {0.0595-in.  diam.  )  orifices.  The  average  total  pro¬ 
pellant  flow  is  2.  45  lbm/sec  at  an  O/F  ratio  of  2.  65.  A  5%  variation  from 
the  average  flow  rate  is  allowed  for  inclusion  in  the  test  results.  The 
chamber  length  can  be  varied  between  8  and  25  inches.  The  chamber  and 
nozzle  throat  diameters  are  2  and  1.  65  in. ,  respectively,  resulting  in  a 
contraction  ratio  of  1.47. 


After  the  start  up  transient  has  been  completed,  and  steady  state 
data  obtained,  the  wave  generator  tube  diaphragm  is  ruptured,  and  a  wave 
propagates  into  the  rocket  motor  thrust  chamber.  The  wave  generator 
tube  is  described  previously  in  this  section. 

Three  high  frequency- response  pressure  transducers  (Fhotocon 
model  no.  352)  are  flush  mounted  in  the  chamber.  They  are  located  2  in. 
from  the  injector  face,  1  in.  from  the  start  of  the  nozzle  contour,  and  mid¬ 
way  between  the  injector  and  nozzle..  Data  are  recorded  on  magnetic  tape 
at  60  ips  and  played  back  at  l  ips  into  a  recording  oscillograph  with  a  paper 
speed  of  25  ips.  Steady  state  pressure  gradients  are  obtained  before  the 
wave  input  by  connecting  twelve  axially  spaced  chamber  pressure  ports  to 
a  single  transducer  through  a  commutating  system. 


The  data  for  the  motors  tested  was  reduced  to  obtain  wave  slope 


histories 


The  wave  slope  is  defined  as  the  ratio  of  wave  pressure 


amplitude  to  wave  rise  time.  The  independent  parameter  selected  for  pre¬ 
sentation  of  the  results  is  the  maximum  value  of  the  steady  state  axial 


pressure  gradient,  (— — )  max.  This  latter  parameter  is  the  maximum 
AX 

slope  on  the  steady  state  pressure  versus  axial  distance  curve  obtained 

from  the  commutator  system.  It  was  chosei.  because  it  readily  describes 

the  type  of  fluid  dynamic  field  through  which  the  wave  propagates.  Varia- 

AP 

tions  in  the  magnitude  of  (  — (max  were  obtained  by  employing  various 

A  X 

combinations  of  injector  configurations  and  chamber  length. 

Figure  4.  14  shows  the  wave  slope  versus  maximum  pressure  gradient 
for  the  incident  wave  at  the  injector  (a),  the  incident  wave  at  the  nozzle  (b), 
and  the  reflected  wave  at  the  injector  (c).  The  shock  driver  tube  pressure 
was  1300  psi  and  the  steady  state  combustion  chamber  pressure  was  approxi¬ 
mately  180  psi  for  all  tests.  Therefore  the  incident  wave  slope  measured 
near  the  injector  face  is  constant  for  most  of  the  pressure  gradients  tested. 

At  the  highest  values  of  pressure  gradient  considered  the  incident  wave  Blope 

AP 

decreases  rapidly.  Since  the  higher  values  of  (— — )  are  produced  closer  to 

AX 

the  injector  face,  wave  slope- gas  dynamic  interactions  have  already  caused 
a  wave  broadening  effect  at  the  first  transducer  location.  The  slope  of  the 
incident  wave  measured  at  the  nozzle  entrance,  (Fig.  4.  14b),  decreases  with 
an  increase  in  the  maximum  slope  of  the  axial  pressure  profile.  Similar 
results  are  obtained  for  the  slope  of  the  nozzle  reflected  wave  measured  near 
the  injector  face;  Fig.  4.  14c. 

More  meaningful  results  can  be  seen  by  normalizing  the  latter  two 
v/ave  slopes  with  respect  to  the  slope  of  the  incident  wave  at  the  injector. 
Figure  4.  15a  indicates  that  the  downstream  propagating  wave  steepens  at 
the  lowest  values  of  the  pressure  gradient  and  broadens  at  the  higher  values. 


i:. 


It  should  be  noted  at  this  point  that  all  maximum  pressure  gradients  les3 
than  3  psi/in.  were  obtained  with  showerhead  injectors  and  all  values 
greater  than  3  psi/in.  were  obtained  with  impinging  injectors.  Normaliza¬ 
tion  of  the  slope  of  the  nozzle  reflected  wave  measured' at  the  injector, 

Fig.  4.  15b,  shows  that  wave  steepening  is  obtained  with  showerhead  injec¬ 
tors  and  wave  broadening  is  obtained  with  impinging  injectors. 

The  ratio  of  the  ordinates  of  Fig.  4.  15b  to  Fig.  4.  15a  were  calculated 
to  determine  the  behavior  of  the  nozzle  reflected  wave  propagating  from  the 
nozzle  to  the  injector.  The  ratio  of  the  wave  slope  at  the  injector  to  the  wave 
slope  at  the  nozzle  is  always  greater  than  unity,  indicating  wave  steepening 
for  all  waves  propagating  unstream  against  an  accelerating  gas  flow. 

The  occurence  of  wavelets  behind  the  initial  wave,  and  the  subsequent 
coalescence  of  these  wavelets  to  produce  a  smooth  fronted  wave  near  the 
nozzle  entrance  has  been  verified  previously,  Fig.  4.  10.  It  is  assumed 
that  these  wavelets  occur  due  to  the  increased  evaporation  rate  and  sub¬ 
sequent  combustion  of  the  fuel  in  an  oxidant  rich  environment.  Some  of  the 
additional  mass-energy  is  used  to  drive  the  wave,  and  part  is  available  to 
increase  the  chamber  base  pressure.  Figure  4.  16  shows  the  ratio  of  base 
pressure  at  the  injector  at  the  arrival  of  the  nozzle  reflected  wave  to  the 
pressure  at  that  point  prior  to  the  incident  wave.  The  data  indicates  that 
pressure  amplification  is  inversely  proportional  to  the  chamber  length  and 
maximum  slope  of  the  steady  state  pressure  curve.  The  amount  of  wave 
initiated  evaporation,  coalescence,  and  subsequent  increase  in  chamber  base 
pressure  is  greater  with  showerhead  injectors  than  impinging  injectors. 
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This  may  be  related  to  the  fuel  drop  size  in  the  chamber.  In  the  impinging 
injector  motors,  the  fuel  drop  size  is  smaller  and  the  conversion  of  liquid 
propellants  to  gaseous  products  is  more  rapid  than  in  the  showerhead 
injector  motors.  Therefore,  there  is  less  mass  and  energy  available  to 
drive  the  waves  or  increase  the  base  pressure  in  the  impinging  injector 
motors. 

For  all  injectors  tested,  the  input  waves  tend  to  be  completely 
attenuated  and  the  hase  pressure  assumes  its  pre-wave  value  after  three 
traverses  of  the  chamber.  At  this  point,  however,  a  second  train  of  longi¬ 
tudinal  waves  occur  in  engines  using  showerhead  injectors.  These  waves  all 
have  rise  times  of  50  microseconds,  independent  of  the  chamber  length.  The 
frequency  of  the  disturbance  (defined  as  the  reciprocal  of  the  time  between 
two  successive  waves  traveling  in  the  same  direction)  corresponds  to  the 
fundamental  longitudinal  Mach  number  compensated  mode.  The  presence 
of  secondary  oscillations  is  not  found  with  impinging  injector  motors. 

In  conclusion,  it  appears  that  wave  steepening  and  pressure  ampli¬ 
fication  are  strongly  coupled  to  the  steady  state  gas  dynamic  flow  field 
through  which  the  wave  must  propagate.  Those  injector- chamber  configura¬ 
tions  which  result  in  rapid  propellant  utilization  and  high  pressure  gradients 
tend  to  inhibit  wave  growth.  Less  rapid  conversion  to  gaseous  products, 
with  an  attendant  low  pressure  gradient,  provides  an  energy  and  mass 
source  to  drive  the  wave  and  amplify  the  base  pressure. 
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SECTION  V 
CONCLUSIONS 

The  purpose  of  this  investigation  in  to  develop  or.  analytical  opray 
comhustion  model  for  a  liquid  bi-propellant  rocket  motor.  Concomitant 
experimental  investigations  of  the  non-oscillatory  steady- state  gas-dynamic 
behavior  and  longitudinal  wave  propagation  phenomena  are  pursued. 

A  two-component  spray  combustion  analysis  is  formulated.  The 

model  considers  separate  evaporation  of  both  fuel  and  oxidizer.  In  addition, 

combustion  gas  properties  are  determined  at  the  local  O/F  ratio  from 

chemical  equilibrium  considerations.  The  analysis  has  been  successfully 

programmed  for  the  I.  B.  M.  7040.  By  removing  the  restrictive  assumption 

4  34 

of  proportional  evaporation  used  in  previous  investigations  ’  a  more 
accurate  description  of  the  aerothermochemical  phenomena  is  obtained. 

The  presently  reported  spray  combustion  model  is  of  immediate 
use  to  the  designer  of  stable  rocket  motor  configurations.  It  can  also  provide 
well  defined  initial  conditions  for  future  analyses  of  combustion  instability. 
The  accuracy  of  the  model  is  dependent  upon  the  accuracy  of  the  assumptions 
used  in  formulating  the  analysis.  Therefore,  more  detailed  studies  of 
injection  processes,  i.  e. ,  jet  fragmentation,  break-up,  droplet  distributions 
and  droplet  shattering  are  recommended  as  a  result  of  the  present  effort. 

The  results  of  the  non-oscillatory  steady  state  experiments  indicate 
that  variations  in  injector  configuration  produce  variations  in  axial  pressure 
history  in  a  constant  area  liquid  propellant  rocket  motor.  These  pressure 
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variations  are  due  to  changes  in  the  mixing  and  shattering  characteristics 
of  different  injectors  and  are  concur  rent  with  unique  evaporation  and  mass 
liberation  profiles.  For  an  injector  that  promotes  rapid  droplet  breakup 
and  intimate  mixing  of  propellants,  the  chamber  pressure  has  a  maximum 
gradient  close  to  the  injector  face  with  a  correspondingly  high  gas  acceleration. 

The  wave  propagation  experiments  subsequently  proved  that  wave 
steepening  and  pressure  amplification  are  strongly  coupled  to  the  steady  state 
gas  dynamic  flow  field  through  which  the  wave  propagates.  Those  injector- 
chamber  configurations  which  result  in  rapid  propellant  utilization  and 
high  pressure  gradients  tend  to  inhibit  wave  growth.  Less  rapid  conversion 
to  gaseous  products,  with  an  attendant  low  pressure  gradient,  provides  an 
energy  and  mass  source  to  drive  the  wave  and  amplify  the  base  pressure. 
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APPENDIX  A 


SPRAY  COMBUSTION  COMPUTER  PROGRAM 

The  program  is  written  for  the  I.  B.  M.  7040  digital  computer,  using 
the  Fortran  IV  language.  The  computational  procedure  is  controlled  by  a 
MAIN  program  which,  in  turn,  controls  a  number  of  subroutines.  Only  the 
MAIN  program  and  those  subroutines  which  form  the  core  of  the  computational 
scheme  are  appended  to  this  discussion.  The  following  is  a  terse  explana¬ 
tion  of  the  entire  program. 

MAIN  PROGRAM 

READ  -  This  routine  reads  in  input  data  for  one  case.  This  data 
includes  the  following. 

CHAMBL  -  length  of  combustion  chamber. 

CONOZL  -  length  of  converging  nozzle. 

FDSHNO  -  fuel  drop  shatter  number;  which  represents  the  number  of 
drops  of  total  equivalent  mass  that  are  used  to  replace  a  single  drop  that  has 
exceeded  the  Weber  number  criteria. 

FJFLOI  -  injected  mass  flow  rate  of  fuel. 

FJVEU  -  fuel  jet  injection  velocity. 

FWTMOL  -  fuel  molecular  weight. 

VARMAX  -  convergence  criteria. 

WEBCR  -  critical  Weber  number. 

NOFDRS  -  number  of  different  size  fuel  drops  used  to  represent 


the  spray. 


NOITRS  -  maximum  ntimber  of  allowable  iterations. 

NOPNTS  -  number  of  points  into  which  the  rocket  motor  is  divided. 

The  next  five  quantities  are  read  in  for  each  of  the  classes  of  fuel 
drops  ranging  from  1=1,  NOFDRS. 

FDCONB(I)  =  concentration  or  percent  of  total  injected  fuel  flow 
represented  by  the  1^  class. 

FDPOSB(I)  =  axial  position  at  which  drops  form. 

FDRADB(I)  =  initial  fuel  drop  radius. 

FDTEMB(I)  =  initial  fuel  drop  temperature. 

FDVELB(I)  =  initial  fuel  drop  velocity. 

All  of  the  above  quantities  which  relate  to  the  boundary  conditions  on 
the  injected  fuel,  (those  containing  F  in  the  symbol),  are  also  read  in  for 
the  injected  oxidizer.  The  READ  routine  also  controls  the  initial  guess  on 
the  values  of  gas  pressure  (PA),  gas  velocity  (VA),  combusted  gas  O/F  ratio 
(EA)  and  gas  temperature  (TA)  at  the  injector.  The  READ  routine  then  calls 
CGGUESS  which  uses  appropriate  functions  to  assume  a  value  for  each  of  the 
last  four  parameters  at  every  point  in  the  chamber-nozzle. 

RDWRIT  -  This  routine  writes  out  the  data  that  was  previously  read 
and  documents  it.  The  initial  combusted  gas  profile  calculated  in  CGGUESS 
is  also  documented. 

DO  400  -  This  statement  controls  the  number  of  iterations  to  be  made 
between  the  liquid  and  gas  systems,, 

COMSET  -  The  last  set  of  combustion  gas  parameters,  pressure, 
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temperature,  velocity,  density,  are  stored  prior  to  the  start  of  the  next 
iteration. 

TITLE  -  This  routine  documents  the  start  of  an  iteration. 

FBSF.T  -  The  method  of  computation  involves  calculating  the  history 
of  a  single  drop  of  each  of  the  drop  classes  and  then  summing  over  all  of 
the  drops  of  all  of  the  classes  to  determine  the  bulk  fuel  parameters,  prior 
to  starting  the  combustion  gas  caculations.  Since  the  bulk  quantities  are 
determined  from  a  summation  procedure,  the  bulk  fuel  properties  at  every 
point  must  be  set  equal  to  zero  prior  to  a  new  calculation  of  droplet  histories. 
FBSET  performs  this  operation. 

DO  200  KLASS  =  1,  NOFDRL  -  This  do  loop  controls  calculation  of 
fuel  droplet  histories  for  each  of  the  classes  of  fuel  drops. 

FDH1ST  -  This  subroutine  calculates  the  history  of  a  single  drop  of 
one  of  the  drop  classes  from  the  point  at  which  droplets  are  formed  to  the 
point  at  which  the  drop  radius  is  zero.  The  actual  solution  of  the  drop 
Eqs.  (9,  10,  11)  is  perfonned  in  FDADV.  The  equations  are  numerically 
integrated  over  a  small  distance  such  that  the  drop  temperature  does  not 
change  by  more  than  10°F.  The  combustion  gas  properties  required  to 
perform  the  calculations  ane  determined  at  each  point  by  interpolation  of 
the  parameters  between  the  bracketing  major  points.  If  the  critical  Weber 
number  is  exceeded,  t<  drop  population,  which  is  initially  unity,  is  multi¬ 
plied  by  the  fuel  drop  shatter  number  (FDSHNO).  In  addition,  the  drop  radius 


is  divided  by  the  cube  root  of  the  shatter  number.  FDADV  also  calculates 
the  energy  of  the  vaporized  fuel  above  the  datum  level,  (FDHES),  as 
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explained  in  Section  II,  Eq.  (26). 

FDWRIT  writes  the  results  of  the  fuel  drop  calculations. 

FMACRO  calculates  the  macro  parameters  for  a  drop  class  by 

multiplying  the  force  on  a  single  drop  (FDFOR),  the  work  done  by  a  single 

drop  (FDFOV),  the  heat  transferred  to  a  single  drop  (FDHET)  and  the  energy 

til 

above  the  datum  level  (FDHES)  by  the  number  of  drops  of  the  I  class  per 
unit  length.  In  addition,  the  flow  rate  of  a  clasB  of  drops  (lbm/sec)  is 
calculated  at  each  point  by  multiplying  the  mass  of  a  single  drop  by  the 
number  of  drops  per  second. 

FMWRIT  writes  the  results  of  the  macro  calculations 
FBSUM  adds  up  the  results  of  the  macro  calculations  from  each 
class  of  drops  at  every  point.  The  following  quantities  are  obtained. 

FBFLO  -  the  total  liquid  fuel  flow  per  second. 

FBFOR  -  the  force  on  all  of  the  liquid  drops  per  unit  length. 

FBFOV  -  the  work  done  on  the  liquid  drops  per  unit  length  per  unit 

time. 

FBIIET  -  the  heat  transferred  to  the  drops  per  unit  length  per  unit 

time. 

FBHEG  -  the  total  energy  of  the  evaporated  fuel  above  the  datum 
level  per  unit  length  per  unit  time. 

FBWRIT  writes  the  results  of  FBSUM  after  all  of  the  fuel  drop 
*  classes  have  been  calculated. 

The  next  nine  statements  in  the  MAIN  PROGRAM,  from  OBSET  to 
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OBWRIT  control  the  calculation  of  oxidizer  histories. 

CGPROP  -  This  routine  controls  solution  of  the  gas  system  Eqs. 

(12,  14,  19,  23).  It  determines  the  ratio  of  evaporated  oxidizer  to  evaporated 
fuel  within  each  interval,  (  JDRAT,  referred  to  as  ELOF  in  Eq.  (25)),  by 
calculating  the  ratio  of  the  decrease  in  OBFLO  to  the  decrease  in  FBFLO. 
The  combusted  gas  O/F  ratio  at  each  point  is  calculated  according  to 
Eq.  (24)  of  Section  II.  The  DO  200  loop  within  CGPROP  calculates  the  com¬ 
busted  gas  flow  at  each  point  by  subtracting  the  bulk  oxidizer  and  fuel  flows 
from  the  injected  flow  rates.  This  DO  loop  also  determines  the  boundary 
condition  on  gas  velocity  by  solving  the  continuity  equation  using  values 
assigned  to  PA,  TA  and  EA  in  CGGUESS  and  the  previously  calculated  gas 
flow  rate.  The  DO  400  loop  controls  calculation  of  gas  properties  from  the 
point  at  which  the  first  drops  form  to  the  nozzle  throat.  The  actual  solution 
of  the  equations  in  performed  in  CGADV. 

CGADV  uses  the  results  of  the  bulk  liquid  calculations  to  numerically 
integrate  the  gas  system  equations.  Values  of  the  energy  coupling  terms, 
i.  e.  ,  BFOR,  BHET,  BHES,  and  BFOV  as  well  as  the  combusted  gas  flow 
(FLO)  are  determined  at  each  minor  subdivision  by  interpolating  between 
the  major  points.  The  number  of  minor  subdivisions  within  each  major 
segment  is  increased  from  8  to  200  in  the  first  four  segments  and  from 
8  to  50  throughout  the  converging  portion  of  the  nozzle.  This  procedure, 
in  the  region  where  large  changes  in  gas  properties  are  probable,  has 
assured  slowly  changing  integrands.  However,  provisions  have  been  made 
in  the  routine  for  inserting  a  test  on  the  magnitude  of  the  change  in  a 


variable,  with  subsequent  step  size  reduction. 

The  Mach  number  of  the  gas  flow  is  calculated  after  convergence 
is  obtained  within  the  minor  subdivision.  If  the  Mach  number  is  greater 
than  one,  the  routine  (CGPROP)  resets  itself,  interpolates  a  new  pressure 
boundary  condition  and  repeats  the  gas  calculation. 

Prior  to  returning  to  the  MAIN  PROGRAM,  CGPROP  will  call 
MAKWON  after  the  DO  400  loop  has  been  completed.  This  latter  subroutine 
calculates  the  Mach  number  at  the  throat  and  compares  it  with  prescribed 
limits.  If  the  throat  Mach  number  is  out  of  range,  a  new  pressure  boundary 
condition  is  interpolated  and  transfer  is  made  to  the  DO  200  loop  within 
CGPROP.  The  range  of  limits  chosen  for  the  throat  Mach  number  is; 

0.  97  <  M  <  1.  03.  The  apparent  3%  deviation  from  the  boundary  condition 
of  M  =  1  is  solely  in  the  interest  of  reducing  computation  time.  Once  a 
pressure  boundary  condition  has  been  found  that  is  consistent  with  the  Mach 
number  at  the  throat  boundary  condition,  transfer  is  made  to  MAIN  PROGRAM. 

CGWRIT  documents  the  results  of  the  last  gas  profile  computation. 

COMPRE  determines  the  maximum  deviation  (VAR)  of  the  gas 
variables  between  successive  iterations. 

The  IF  (N9.  LT.  0)  statement  forces  a  return  to  the  drop  calculations 
if  the  pressure  boundary  condition  has  been  changed  since  the  last  drop 
calculations. 

The  IF  (VARMAX- VAR)  statement  tests  for  convergence  to  a  solution. 

If  convergence  has  not  been  attained,  the  drop  calculations  are  entered  using 
information  obtained  from  the  last  gas  profile  calculations.  When  the 
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convergence  criteria  is  satisfied,  the  program  leaves  the  major  DO  400 
loop. 

MAKWON  -  This  routine  lias  been  described  previously.  Calling  it 
at  this  point  in  the  program  is  actually  redundant  since  the  Mach  number 
at  the  throat  has  previously  been  tested  prior  to  exit  from  CGPROP. 

GO  TO  100  -  This  statement  permits  the  calculations  to  be  repeated 
if  more  than  one  set  of  droplet  boundary  conditions  were  initially  specified. 

The  following  list  of  function  subroutines  are  required. 

CHAREA  -  cross-sectional  area  as  a  function  of  axial  distance. 

FDENS  -  fuel  density  as  a  function  of  temperature. 

FHEVAP  -  fuel  heat  of  vaporization  as  a  function  of  temperature. 

FLSPH  -  liquid  fuel  specific  heat  as  a  function  of  temperature. 

FSTEN  -  fuel  surface  tension  as  a  function  of  temperature. 

2 

FVPRE  -  fuel  vapor  pressure  (lbf/ft  )  as  a  function  of  temperature. 

FVSPH  -  fuel  vapor  specific  heat  as  a  function  of  temperature. 

A  similar  set  of  property  subroutines  must  be  specified  for  the  oxidizer. 
These  are  titled  and  called  in  the  program  by  replacing  F  with  O.  The  data 
may  be  given  as  an  equation  or  it  may  be  specified  in  tabular  form,  with  the 
value  of  the  property  determined  by  an  interpolation  subroutine  (TERP).  An 
example  of  the  former  case  (OVPRE),  and  the  latter  (FHEVAP)  are  given 
in  the  sample  listings. 

The  diffusivity  D^.,  defined  in  Section  II,  is  given  for  the  fuel 
(FDIFUS)  and  the  oxidizer  (ODIFUS)  as  a  function  subroutine.  The  arguments 
arc  gas  temperature  and  pressure.  Diffusivity  was  included  as  a  separate 


subroutine,  rather  than  in  the  FDADV  routine,  since  it  is  a  function  of  the 
particular  propellants  being  analyzed. 

CGWMOL  -  combustion  gas  molecular  weight  as  a  function  of  tem¬ 
perature  and  O/F  ratio 

CGSPHT  -  combustion  gas  specific  heat  as  a  function  of  temperature 
and  O/F  ratio 

CGVISC  -  combustion  gas  viscosity  as  a  function  of  temperature  and 
O/F  ratio 

CGENTS  -  heat  of  reaction  of  the  propellant  combination  per  pound  of 
products,  measured  at  the  base  temperature  as  a  function  of  O/F  ratio. 

In  addition  to  the  above  data,  information  must  be  available  concerning 
the  enthalpy,  temperature  and  O/F  ratio  of  the  combusted  gas.  The  base  level 
for  enthalpy  must  of  course  be  consistent  with  that  used  in  CGENTS,  In  the 
present  analysis  enthalpy-temperature  relations  are  tabulated  for  various 
O/F  ratios  ranging  from  zero  to  100.  A  double  interpolation  routine  (DP) 
determines  either  enthalpy  or  temperature  whenever  the  O/F  ratio  and  either 
of  the  preceding  variables  is  specified.  The  arguments  associated  with  DP 
are  O/F  ratio,  enthalpy  and  temperature,  in  that  order.  If  for  example, 

T  and  O/F  are  known,  and  enthalpy  is  to  be  determined,  the  following 
instructions  are  required,  (see  CGADV  listing). 

HH  =  0 

CALL  DP  (OF,  HH,  T) 

The  enthalpy  will  be  interpolated  and  stored  in  HH.  In  general,  the 
argument  to  be  interpolated,  is  set  equal  to  zero  prior  to  calling  DP. 
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0  SIBFTC  BOB  DECK 

1  DIMENSION  CG06NI  129) ,CCFLO( 129  )  , CGPOS ( 129 ) , CGPRE ( 1 29 > ,CGTEK(129 )  , 
i  CGVEL (  129)  jCGVkMCI  129) 

2  DIMENSION  F OR  A T I  1 29 ) , E I R A  1(1 29 ) 

3  DIMENSION  FBFLOI  129 ) , FBFOR 1 129  > , FBFG V ( 1 29 ) , F BHFS ( 1 29 ) , FBHET ( 129 ) , 

1  FBPUPI  129)  ,F8VEL ( 129) , 

2  FDFoR l 129)  ,FDFCV< 129) , F  CHE  S (  1 29  ) , FCHE I II 29 ) , FDPOP ( 129 )  , 

3  F OPUS (  129)  ,  FDR AC  I  129) ,FCTEM(  1  29  )  , FDT I M 1 1 29 >  » 7 DVEL 1129), 

4  FMFLOI  129  )  ,FMFCR( 129) , FMFOVt  129)  ,FMHtS( 129 ) , FMHE T ( 1 29 ) , 

5  FMPOP ( 129 ) 

4  DIMENSION  FUCONB ( 33 ) , FDPGSB (  33  ) , FDR ADB ( 33  )  .FCTt'MB (  33  ) , FDVELB { 33 ) 

5  DIMENSION  CBFLOI 129) ,CBFOR( 129 > , OBFOV ( 1 29 ) , OBHES ( 129 ) »OBHET ( 129 > , 

1  OBPOPI  129)  ,OBVEL( 129) , 

2  ODFURI 129)  ,OOFCV( 129) , ODHE S 1 1 29 ) , CDHET ( 129 ) ,ODPOP ( 129 ) , 

3  OOPUSI  129)  ,ODRAC(  129)  , OD T EM  ( 1 29  » ,C!CT  IM  1129  )  ,ODVEL  (129), 

4  OMFLOl 129) .OMFCRI 129) .OMFOVI  1 29 ) , OHHES 1 1 29 ) , OMHE TC 129 ) , 

5  OMPOP 1129) 

6  DIMENSION  CDCONB (33 ) , ODPCS B ( 3 3 ) , OCR ADB ( 33 ) ,OCT EMU ( 33 ) pODVELB ( 33 ) 

7  DIMENSION  PAGE  I  I ( 12 ) , PROGT I ( 12 ) 

10  DIMENSION  SCGQEM 129) ,SCGFLO( 129)  ,SCGPRE( 129),SCGTEM(129), 

1  SCGVEL (  129  ) 

11  COMMON  CGDEN,CGFLG,CGPCS,CGPRE»CGTEM,CGVEL,CGhMO 

12  COMMON  EDRA1  ,E  IRAT 

13  CGMMON  FBFLO, FBFOR, FBFO V , F BHES , FBHET , FBPCP .FBVEL , 

1  FOFOK,FDFCV,FCHCS,FOHET,FOPOP,FDPOS,FCRAD,FDTEMtFDTIM,FDVEL, 

2  FMFLO,FMFCR,FMFCV,FMHES,FMHET ,FMPCP 

14  COMMON  FDCGNB,FOPQSB,FCHADB,FCrEMB»FCVELB 

15  COMMON  OBFLO, OBFCR.CBFCV, OBHES, OBHlT.OBPCP.OBVEL, 

1  ODFGR,GDFCV,OCHES,OOHtT,OCPCP,ODPOS,ODRAD,ODTEM,ODTIM,ODVEL, 

2  OMFLO,OMFCR,OMFCV,CMHES,OMHET,OMPCP 

16  COMMON  0DCCNB,GCP0SB,CCRAD6,CCTEMB,0CVELB 

17  COMMON  PAGE  T I , PROG  1 1 

20  COMMON  SCGDEN,SCGFL0,SCGPRE,SCG1EM,SCGVEL 

21  COMMON  CHAMBL .CONQZl , GINOZL,FOSHNO,FJFLOI,FJVELI»FHTMOL»GENTHS» 

1  GHrMOL.ODSHNO.OJFLCI ,CJVELI ,CWTMOL,VARMAX,hEBRCR 

22  COMMON  NOFDPS,NCITRS,NGQDRS,NCOPTS 

23  COMMON  LT IN.LTOUT 

24  EQUIVALENCE  ( NOFP IS , NCGPTS , NCOPTS , NOPNTS > 

25  COMMON/CLASS/KLASS 

26  C0MMCN/FCRCE/N9 

27  LT I N=5 

30  LT0UT=6 

31  100  CONTINUE 

32  CALL  REAO 

33  CALL  RDWRIT 

34  CUMMON/BOUDRY/  P A , T A , V A , ELM , E A 

35  COMMON/ TEST/  NA ,NB , PM  I N , PMAX 

36  OC  40C  I T= 1 ,NOI TRS 

37  IFIIT.GT.2)  GO  TO  IG1 

42  PMIN=PA 

43  PMAX=PA 

44  NA=0 

45  NB=C 

46  101  CONTINUE 

47  CALL  COMSET 
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50 

CALL  TITLE  1 1 T 1 

51 

CALL  FOSET 

52 

00  200  KLASS= 1 « NOFCRS 

53 

CALL  FUH 1ST (KLASS) 

54 

CALL  FOUR 1 1 ( KLASS ) 

55 

CALL  FMACRO (KLASS ) 

56 

CALL  FMHRITCKLASS) 

57 

CALL  FBSUH ( KLASS i 

60 

200 

CONTINUE 

62 

CALL  FBWRIF 

63 

CALL  CBSET 

64 

DO  300  KLASS  =1,N0CDRS 

65 

CALL  ODH 1ST (KLASS ) 

66 

CALL  QDWR I T ( KLASS ) 

67 

CALL  OMACRO (KLASS ) 

70 

CALL  OH WRIT (KLASS ) 

71 

CALL  OBSUMIKLASS) 

72 

300 

CONTINUE 

74 

CALL  0BWR1T 

75 

20002 

CALL  CGPROP 

76 

CALL  CGWRIT 

77 

CALL  COMPRECVARJ 

100 

IFCN9.LT. C)G0r04GC 

103 

I F ( VARMAX-VAR )  400,400,500 

104 

400 

N9=l 

106 

500 

CONTINUE 

107 

CALL  MAKWONCCGVELCNCPNTSI , CGTEM ( NOPNTS ) ,N9 1 

no 

IFCN9.LT. OIGOT020G02 

113 

GO  TO  100 

114 

END 

1IBFTC  FJHIST  D^CK 

SUoRQijT  I  Nfc  FfJM  [Ml  K.LA'jh  ) 

COMMON / F  H  f  A  T  /  S  H  I  ‘ !  T 

DIMENSION  CGDhMI 1 29 ) , CO F LO ( 1 29 ) ,CGPOS< 129)  ,CGPRE (  129  )  ,  C  G  T  E  M  ( 129)  » 
1  CGVELI  129)  ,CGWMO( 120 ) 

DIMENSION  FORATI 120 ) ,f  f pat ( 1 2Q) 

DIMENSION  Fi'FLO  I  129  )  ,  FoFOR  (  1  29  )  »FdFCV<  129)  *  F uHE S  !  129)  ,FBHET(  129)  » 

1  FRPOP{  129),F13VEL(  120)  , 

2  FDFUR  (  129  )  ,  FDFOVI  129  )  ,FDhES(  129  )  >  F  OiHET  (12  9)  .FDPOP  ( 129  )  , 

3  FDPOSI  129), FUKADI 129 )  ,FUTtM( 129 ) ,FDT I  MI  129 )  ♦  r  D  V  E  L  <  129  )  , 

4  FMFlOI 129) , FMFOR  (  129) ,FKFOV< 129) , FMHES (  129) , FMHET (129)  , 
s  F’^POP  (120) 

DIMENSION  FDCONB ( 33 )  ,FDPOSB(33) ,FDRADB ( 33 ) *  FDTEMB (33)  ,FDVELB(33) 
DIMENSION  OBFLO (129)  ,OBFOR( 129) .OBFOVI129) ,OUHES( 129)  ,0BHET(129) * 

1  OBPOP  (  129)  ,0l!VEL  (  129)  , 

2  ODFOR ( 129 ) ,OOFOV<  129 ) ,ODHES(  129 ) ,UOHlT  (  12  9 ) ,ODPOP( 129  )  » 

3  ODPOS ( 129 ) ,ODRAD(  129) ,ODTEM( 129 ) .OUT  INI  129 ) ,ODVEL( 129  )  , 

4  OMFLOI 129 ) , CM FOR  (  129 )  ,OMFOV( 129 ) ,OMHES(  129 ) .OMHET ( 129)  , 

5  OMPOP (129) 

DIMENSION  ODCONB ( 33 )  ,ODPOSB(33) *ODRADb ( 33)  ,0DTEMB(3?) ,ODVELB(33) 
DIMENSION  PAGETI ( 12) ,PROGTI ( 12) 

DIMENSION  SCGDENI 129 ) »SCGFLO( 129 ) *SCGPRE(  129 ) , SCOT  EM ( 129) * 

1  SCGVEL ( 129 ) 

COMMON  CGDEN  *  CGF  L0»  CGPOS  .CGPRt  »CoTEM*CGVt.L  .CGrtMU 
COMMON  EDR AT  *  E I  RAT 

COMMON  FBFLO*FbFOR«FoFOV  »  FbHES  *  FbFiE  T  ,  FBPOP  *  FoVEL  » 

1  FDFOR ,  FDFOV. FDHESt  FDHET ,  FDPOP ,  FDPOS » FDR  AD  *  FDTEM  *  FDT  IM»  FDVEL  ,  ‘ 

2  FMFLO, FMFOR, FMFOV,  FMHE S, FMHE T .FMPOP 
COMMON  fdcONB, FDPOSR,  FDR ADB, FDTEMB, FDVELB 
COMMON  OB  FLO, OBFOR ,OBFOV ,OBHE S , OBHET , OBPOP,  OBVEL, 

1  ODFOR, ODFOV,OUHES,ODHET ,ODPOP .ODPOS , ODR AD , ODT EM ,ODT I M .ODVEL , 

2  OMFLO.OMFOR .OMFOV.OMHES. OMHET. OMPOP 
COMMON  ODCONB*  ODPOS0 ♦ ODR ADB ♦ ODT EMB .ODVELB 
COMMON  PAGETI ,PROGT  I 

COMMON  SCGDEN,SCGFLO,SCGPRE,  SCGTEM  .SCGVl.L 

COMMON  CH AM3L  * CONOZL  ,  D  I  NOZ L  » FDSFINO ,  F  JF'LOI  ,F  JVELI  ,  FWTMOL  .GENTHS , 

1  GWTMOL.ODSHNO.OJFLOI  .OJVELI  ,  OWTMOL  ,VAR‘*AX  ,  WEBRCR 
COMMON  NOFDRS.NO I TRS , NOODRS , NOCPT S 
COMMON  LTIN.LTOUT 

EQUIVALENCE  <NOFPTS,NOGPTS.NOOPTb.NOPNTS> 

IFJL0=1 

IFJHI=1 

DO  300  1  =  1  ,NOPNTS 

IF  (FDPOSI  I  )  -FDPOSB  (  KLA  SS  )  )  2  00 , 3lu  ,  3Ui< 

200  CONTINUE 
I FUH  I  =  I 
3uO  CONTINUE 

IFDLO= IFJHI+1 
IFDHI=N0PNTS 
DO  400  I  =  I  F JLO , I F Jri  I 
FDFOR ( I ) =0. 

FDFOVI I ) =0. 

FDHES ( I)=0. 

FDHET ( I ) =0. 

FDPOP ( I ) =0. 
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FDRADI I ) =0  , 

FDTEMI I ) =0. 

FDTJMI  I  )=FDPOS(  I  1/FJVELI 
FDVELI I ) =0  o 
400  CONTINUE 
P0P=1  o 

P0S=FDP0S ( I FJHI ) 

RAD=FDRADB(KLASS> 

TEM=FOTEMb(<LASS) 

T IM=FDT I M ( I F JH I ) + ( P0S-FUP0S < I FJHI ) I/FJVELI 

VEL=FDVELB(KLASS1 

SH I  NT-0 . 

DO  800  1  =  I FDL0»  IFDHI 
I F ( RAD  I  5  00  »  500  1 600 
500  CONTINUE 

FDFOR ( I ) =0. 

FDFOV ( I ) =0. 

FDHES ( I ) =0. 

FDHET ( I ) =0. 

PBP0D ( I ) =0  , 

FDRADI I ) =0. 

FDTEMI I ) =0. 

FDTIMI I ) =0 . 

FDVELI I ) =0 • 

GO  TO  800 
600  CONTINUE 

POS  =  FDPOS I  I  — 1 1 

CALL  FDADVI I , POP 9POS 9 RAD 9TEM » T I M 9 VEL i 
800  CONTINUE 
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subRour  i  wf  rrvnv  <  nxtpt  ,pop,pos*kad»tem,tim,vel) 
co^mon/cl ass /k  lass 
COMMON/ FHF AT/  5  H  J  N  T 

DIMENSION  COD IN (  I  29 )  . C OF LO ( 1 29 ) tCOPOS ( 1 29 ) tCGPRt ( 129 ) .CGTEM ( 129 )  f 
1  CGVEL (129) »  CGftMO (129) 

DIMENSION  EDRAT ( 129 )  , El  RAT (129) 

DIMENSION  FOFLOI  129)  .EBFORI  129)  .FbFOVt  129)  ,FBHES(  129)  ,FBHET(  129)  * 

1  FBPOP (1 29 ) , FBVEL ( 129 ) , 

2  FDFOR ( 129) , FDFCV ( 129) *  FDHES ( 129) «  FDHET (  129) , FDPGP 1129) , 

■>  FDPCS  (  129)  ,  FDR  ADI  129)  ,FDTEM( 129)  ,FDTIM( 129 J  .FDVELU29)  , 

4  FMFlOI  129)  , FKFOR (  129)  ,FMFOV<129)  .FiMHCSI  129)  ,FMHET(129)  . 

5  FMPOP (129! 

DIMENSION  FDCOiMbl 33 ) .FUPoSol 33! .FDRADol js)  » FDTEKo ( 33 ) *FDVELb  (  33 ) 
DIMENSION  03FL0 ( 129 ) .OBFOR (129) .OBFOV (129) ,OBHES( 129) .OBHETI  129) * 

1  OBPOPt 129) ,CBVEL< 129) , 

2  ODFOR  (  129  )  ,ODFCV(  129!  .ODFiESI  129)  *OuFiET(  129)  ,ODPOP(  129  )  , 

3  ODPOS ( 1 29 ) »ODRAD ( 129)  ,ODTEM( 129 ) *CDT INK  129)  ,ODVEL(129), 

4  OMFLOI 129 )  , OMFOR ( 129 ) .OMFOVI 129) .CMHESI 129 ) ,OMHET (129)  * 

9  OMPOP (129) 

DIMENSION  ODCONb ( 33 ) .ODPOSB(33) .OURADIK  33) * ODTEMB ( 33 ) -ODVELb(33) 
DIMENSION'  PAGLT 1(12) .PROGT  1(12) 

DIMENSION  SCGDEiN  (  129)  ,SCGFLO(129)  *SCOPRE(  129)  , SCGTEM ( 1 29 ) ♦ 

1  SCGVEL (  129 ) 

COMMON  CGDEN.CGFLO.CGPOS.CGPRE, CGTEM,  CGVEL, CGwMO 
COMMON  FDR  AT .FIR AT 

COMMON  FBFLO.FBFOR.FBFOV.FBHES.FBHET. FBPOP. FBVEL. 

1  FDFOR , FDFOV. FDHES  *  F DHtT ,FDROP»FDPOb»FDRAD  *FDTEM»FDT IMtFDVEL ♦ 

2  FMFLO.FMFOR .FMFOV.FMHES.FMHET, FMPOP 
COMMON  FDCONB. FDPOSU , FDRADB  ,  FDTEMo , FuVELB 
COMMON  OBFLO.OPFCR.OBFOV.OBHES.ObHET.OBPOP.OBVEL. 

1  ODFOR, ODFOV,ODHES»ODHET  .ODPOP .ODPOS » ODRAD .ODTEM.ODt I M.ODVEL , 

2  OMFLO.OMFOR .OMFOV. OMHES  »OMHET .OMPOP 
COMMON  ODCONB, ODPOSB , ODRADB , ODTEMB, ODVELB 
COMMON  PAGET  I »  PROGT I 

COMMON  SCGDEN.SCGFLO.SCGPRE,  SCGTEM, SCGVEL 

COMMON  CHAMBL.CONOZL.DINOZL.FDSHNO.FJFLOI  ,F  JVELI  .F'aTMOL.GENTHS, 

1  GWTMOL .ODSHNO.OJFLUl .OJVELI ,0*TK0L . VAKMAX » rtEBRCR 
COMMON  NOFORS.NO I TR S , NOODRS , NOOPTS 
COMMON  LT IN.LTOUT 

EQ'J  I  VALENCE  (  NOFPTS  , NOGPTS  .NOOPTS  , NOPNT  S  ) 

GRCON=32 . 2 
UGCON= 1 54  5 • 

NSUBDV=8 
NN  =  1 

126  DO  1  GOG  J  =  NN  »  NSUBDV 

DLPOS=(FDPOS(NXTPT ) -FDPOS ( NXTPT-1 ) ) /FLOAT ( NSUBDV ) 

IF (RAD.LT.FDRADcl KLA5S) / < 10. *POP**. 33333 )  )  GO  TO  700 
FRACT=(POS-FDPOS( NXTPT-1 ) ) /( FDPOS ( NXTPT ) -FDPOS ( NXTPT-1 ) ) 
GEIR=EIRAT (NXTPT-1 ) +FRACT* ( FIR AT ( NXTPT ) -E IRAT ( NXTPT-1 ) ) 
GPRE=CGPRF(NXTPT-1 )+FR ACT* (CGPRE( NXTPT )-CGPRE( NXTPT-1 1 ) 

GT EM =C GTEM ( NXTPT-1 ) +FRACT *(CGTEM( NXTPT )-CGT EM (NXTPT-1 ) ) 
GVEL=CGVEL(NXTPT-1 )+FR ACT* (CGVEL (NXTPT )~CuVEL( NXTPT-1 ) ) 

GWTM  =  CGWM0(NXTPT-1  ) +FRACT*  (  CGWMOI  NXTPT  ) -CGWf-  0  (  NXTPT-1  )  ) 

TEMBAR= ( TEM+GTEM ) /2 • 
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GSPH=CGSPHT< TEMBAR, GE1R ) 

GV  I S=CGV I  SC  t TEME  *R»GWTM ) 

REYNUM=2 .  »RADftABS( GVEL-VEL ) *GPRE*GWTM 
X  /(UGC0N«TEMBAR*GVIS) 

SCHNUM  =  GVI  5/  (FDIFUSI TEMBAR ,GPRE  )  *GPRE*GWTM 
X  / (UGCON& TEMBAR)  ) 

COMAST=  <  2e+u.6«SCHNUM**  (  1 • /3.  ) *REYNUM** { 1 . /2 .  ) 

X  *FDI  FUS  (  TEMBAR, GPRE  )  »FWTMOl/  (  UC-C0N*2.*R AD# TEMBAR  ) 

DARFA=4„*3.1415926»RAD'*2 
DLTIM=DLPOS/VEL 
XDDMDT  =DDMDT 
IF <  TEM.EQ.O.  )•  GO  TO  730 
IF(GPRE.LT.FVPREITEM)  >  GO  TO  90 
ODMDT  =  DAREA*COMAST*GPRE*ALOG  ( ABS  ( GPRE  /  ( GPRE-FVPRE  ( TEM )  )  )  ) 

0MEGA  =  AMAX1  (DDMDT, XDDMDT  J  . . 

GO  TO  98 
90  CONTINUE 

LOGICAL  DOING 
BOING=DLTEM.NE.O. 

IFlBOlNGoAND. NSUBDV.LT. 200)  GO  TO  99 
DDMDT  =  OMEGA 
GO  TO  600 
98  'CONTINUE 

H=J2.+U.525*SQRT(REYNUM)  )  *GSPH*GV  I  S/ ( 1 . 34*RAD ) 

2  =DDMDT»FVSPH( TEMBAR ) /  (  H*DAREA ) 

I  F ( 2-1 oE— 4  1  200,300,300 
200  CONTINUE 
2=1. E-4 
GO  TO  500 
300  CONTINUE 

I F ( 25 « -Z  )  400,530,500 
400  CONTINUE 

~ .  T=25. - 

500  CONTINUE 

HDNET  =  H*DAREA*( GTEM-TEM i  *Z/ ( EXP ( Z ) -1 o  ) 

HD  =  HDNET -DDMDT *FHEVAP  (  TEM ) 

DVOL=4./3.*3.1415926«RAD**3 
HCAPD=FDENS  <  TEM )*DVOL*FLSPH ITEM) 

DLTEM- tHO/HCAPD) *DLT  IM 
I F ( ABS ( DLTEM ) .GE.10. )  GO  TO  100 
IF (TEM+DLTEM)  700,700,601 
600  DLTEM=0« 

601  CONTINUE 

CODRAG=27./REYNUM««(0.84) 

ACCEL=  C  3  r,  /  8  .  )*C0DRAG»  ( GVEL-VEL ) *AbS ( GVEL-VEL  ) 

X  * (G*TM»GPRE/ (UGCON* TEMBAR )  )  /  ( FDENS ( TEM) *RAu ) 
DLVEL=ACCEL*DLTIM 
I F ( VEL+DLVEL )  73.;  ,730 ,602 

602  CONTINUF 

DLRAD=-DDMDT*DLTIM/(DAREA*FDENS( TEM) ) 

X  -RAD*  (  FDENS  (  TEM+DLTEMT-FDENS  (  TEM)  )  / 1 3  .^FDENSlTEMTI 

I F ( RAD+DLRAD  )  700,700,603 

603  CONTINUE 
GO  TO  750 

700  CONTINUE 
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-0 . 

X  » i  *  '■'  D  T  =■  « 

PCPO. 

RA0=O. 

DLRAD=„  . 

T  Li’i=  -  • 

DLTEM- J  . 

TIK=j. 

OLT I M  =  0 . 

VEL  =  0 . 

DLVEL=0. 

75G  CONTINUE 
90C  CONTINUE 
XHEP=HEP 
XOHINT=SHINT 
XPOS=POS 
XPOP=PCP 
XRAD=R AD 
XTEM=TEM 
XT  I M=T I M 
XVEL=VEL 

XX=FDENS (  TEM ) *RAD**3 
POS=POS+DLPOS 
RAD  =  R AD  +  DL  RAD 
TEM=TFM+DLTEM 
TIM=T  IM+rLl'IM 
VEL=VEL+DLVEL 
YY=FDENS(TEM)*RAD  #*3 
SHINT=SHINT+FLSPH(TEM)*  DLTEM 
dt=dltim 

HEP  =  HEP + ( SH I NT  +  FHEVAP ITEM)  1*1.33333*3.  1 4  1  59*ABS ( XX-YY ) *POP/DT 
1  *l./FLOAT(NSUBDV) 

WBNUM=GWTK*GPRE* (GVEL-VEL  )**2*2.*(RAD  ) 

X  /  ( UGCON*TEvIBAR*FSTEN  (  TEM  1  1 
IF(WEBRCR-WBNUM)  800.800,1000 
8C0  CONTINUE 

RAD=RAD/FDSHNO**, 3333333 

pop=pop*fdshno 

OMEGA=RAD**2/XRAD**2*OMEGA 
DDMDT  =OMEGA 

1000  CONTINUE 

1001  CONTINUE 
FDHESt  NXTPT ) =HEP 
HEP  =  0 

FDFOR  (NXTPT 1 =DVOL*FDENS ( TEM ) *DLVEL / ( ULT I M*GRCON 1 

FDFQV (NXTPT 1 =FDFCR ( NXTPT  5*VEL 

FDHET ( NXTPT ) =HDNET 

FDPOP ( NXTPT ) =POP 

FDRADt NXTPT )=RAD 

FDTEM(NXTPT) =T  EM 

FDTIM(NXTPT) =TIM 

FDVEL(NXTPT)=VEL 

RETURN 

99  IF ( J.EQ.l  )TEM=TEM-AoG(ULT£M) 

100  NN=J*2-3 
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0LP0S=DLP0S/2. 

NSUnnV=NSUBDV*2 
rF(NN.LT.P)  G 0  TO  113 

pos=xpos 

POP=XPCP 
RAD=XRAO 
TEM=XTEM 
T  !  M  =  X  T  I M 
VEL  =  XVFL 
HEP  =  XHFP 
SH I NT  =  XSH INT 
GO  TO  126 
113  NN= 1 

GO  TO  126 
END 

$  I  5FTC  ODADV  D6CK 

SUBROUT l NEODADVI NXTPT .POP .PCS, RAD, TEM,  TIM, VEL) 

COMMON/CLASS/KLASS 

COMMON/FHFAT/  SHINT 

DIMENSION  CGDEN  (  129 ) . CGFLO (  129) .CGPOSI 129)  .CGPREl 129 ) ,CGTEM(  129 ) 
1  CGVELI129).  CGWMO (12  9) 

DIMENSION  EDRAT ( 129 ) . L I  RATI  129) 

DIMENSION  FBFLC ! 129)  ,FBFOR( 1 2 o ) .FBFOVI 129)  , FbHESI  129) .FBHETI  129) 

1  FBPOPl 1  29  )»  FBVEL ( 129)  . 

2  FDFORl  129)  .FDFCVI  12°)  .FDHESU29)  » FDHET  (  129)  .FDPOPI129), 

3  FDPOSI  129)  ,  FDR  AD  I  1  29  )  .FDTEMI 129  )  .FDT  I  MI  129  )  .FDVELU29)  . 

A  F.MFLOI  129  )  ,  FMFOR  I  129  )  .FMFOVI  129)  .  FMHES  I  129)  .FMHETI  129)  , 

5  FMPOP (129) 

DIMENSION  FUCONEI 33 ) , FUPOSbl 33) .FDRAUbI  33)  , FuTEMbl 33) ,F DVELdI  33 ) 
DIMENSION  OBFLOI 129)  ,UBFOR(129! ,ObFOV(129) .OBHESI 129) .ObHETI  129) 

1  OBPOP(129),OOVEL: 129) , 

2  ODFORI 129) .ODFOVI  129)  .ODHFS I  1 29 ) , ODHET I  129) .ODPOPI129)  , 

3  ODPOSI 129) .ODRADI 129) .ODTEMI 1 29 ) . ODT I M C 1 29 )  .ODVEU129) , 

4  OMFLOI 129)  . OMFQR (129) .OMFOVI 129) .OMHESI 129 ) .OMHET 1129). 

5  OMPOP (129) 

DIMENSION  ODCONB (33) , JoPOSB(33) .ODRADb ( 33)  .ODTEMB ( 33 ) ,ODVELB(  33 ) 
DIMENSION  PAGETI ( 12) .PROGTI ( 12) 

DIMENSION  5CGDEN ( 129  )  .SCGFLO( 129 ) .SCGPRE ( 129 ) .SCGTEMI 129)  . 

1  SCGVEL (  129 ) 

COMMON  CGDEN, C GF LO, CGP OS, CGPRE.CGTEM.CGVEL. CGWMO 
COMMON  FDRAT.FI RAT 

COMMON  FR FLO , FG FOR , FBF  OV , FBHES ♦ FBHET  » FBPOP , FBVEL , 

1  FDFOR.FDFOV.FDHES.FOHET.FDPOP.FDPOS.FDRAD.FDTEM.FDTIM.FDVEL , 

2  FMFlO  ,  FMFOR  ,  F F 0  V  ,  FMHES  ,FMHE T  .FMPOP 
lOMMON  FDCONB »  FDPOSrf , FDRADB , FDTEMb , F jVELB 
COMMON  OB  FLO  »ObF  OF;  »  DoF  GV  .ObrIE  S  »Ot>Ht  T  »  UbPUP  .UUVEL  , 

1  ODFOR  ,  OOFOV  ,0  JHE  S  »  C’DHE  T  .UUPOP  .OUPOS  »  ODR  AD , 0D  T  tM  » OUT  I M , ODVEL  , 

2  omflo.cmfok  .c.'FD'V  ,'omhes  .qmhet  .ompgp 

COMMON  ODCONB,  RDP05P  .ODRADB  .ODTEMfc .ODVELL* 

COMMON  PAGET  I .PROGT I 

COMMON  SCGDEN.SCGFLO, SCGPRE, SCGTEM, SCGVEL 

COMMON  CHAMBL , CONOZL , D I NOZL , FDSHNO , F JFLO I , F JVEL  I  .FWTMOL .GENTHS , 

1  GWTMOL.ODSHNO.OJFLOI , OJVEL I , OtoTMOL .VARMAX.WEoRCR 
COMMON  NOFDRS.NCI TRS.NOODRS.NOOPTS 
COMMON  LTIN.LTOUT 


r  0  U  I  V  A  L  r  N  c  F  ( r  i  0  F  °  T  G , : .  ^  G  p  t  S  ,  N  0  0  p  T  G  ,  N  OP  NT  S ) 

GRCON=32.2 
'JGCON  =  1545. 

NSUBDV  =  2 
N  N  *  1 

126  DO  lOuO  J  =  NN , KSU3DV 

DLPOS= ( FDPOSt  NXTPT )- F DPC3 ( NXTPT - 1 ) ) /FLOAT ( NSUBDV ) 
IF(RAD.LT.ODRAD:5(KLAGS  )/(  10.  *POP**.  33333)  )  GO  TO  700 
FRACT  =  (  POS-ODPOS  ( TiX T PT  - 1  )  )  /  (CDPOSC  NXTPT)  -ODPOS  (  NXTPT -1  )  ) 
GE I R  =  c I RA  T  fiXTPT-1 )4TRACT*(EIRAT I NXTPT ) -E I  RAT ( NXTPT -1 ) ) 
GPRr  =  rGn>^F(  NXTPT- 1  )  +  FR  ACT*  (  CGPRF  (  NXTPT  )  -CGPRE  (  NXTPT- 1  )  ! 
GTCM=CGTEM ( NXTPT- 1 ) +  FR ACT* ( CGTEM ( NXTPT 1 -CGTEM < NXTPT -1 ) ) 
GVEL=CGVEL (NXTPT- 1 ) +FR ACT* ( CGVEL ( NXTPT ) -CGVEL ( NXTPT -1 ) ) 
GWTM=CGWMO( NXTPT- 1 ) +FR ACT* ( CGWMO ( NXTPT ) -CGWMO ( NX TPT-1 ) ) 

T  EMbAR= ( TEM+GTEM) /2. 

GSPH=CGSPHT ( TEMBAR ,GE I R ) 

GV I S =CGV I  SC ( T EM BAR , GW TV ) 

REYNUM=2. *RAD*ASS(GVEL-VEL ) *GPRE*GWTM 
X  / <UGCON*TEMBAR»GVIS) 

SCHNUM=GV  1  5/  (03  I  FUG  (  TlMliAR  ,GPRE  )  *GPRE*GWTM 
X  /  (  UGCON*T  EMtiAR  )  ) 

COMAST = ( 2 .+J.6*SCHNUM»* { 1 . / 3  « )*REYNUM** ( x./2  .  )  ) 

X  *  0  D I F  U  S ( TEMBAR.GPRt ) *OWTMOL/< UGC0N*2.*RAb*TEMBAR ) 

DAREA=4.*3.1415926*RAD**2 
XDDMDT  =DOMDT 
1 F ( TEM.EO.O. )  GO  TO  700 
IF(GPRE.LT.OVPRE(TFM)  )  GO  TO  90 
DDMDT =DAREA*C0MA3T*GPRC*AL0G( ABS ( GPRL/ ( GPRE-OVPRE ( T EM ) ) ) ) 
0MEGA=AMAX1 (DDMDT. XDDMDT ) 

GO  TO  98 
90  CONTINUE 

LOGICAL  BOING 
BOING=DLTEM.NE.O. 

IF(  BOING. AND. NSU5DV.LT. 200)  GO  TO  99 
DDMDT  =OMEGA 
GO  TO  600 
98  CONTINUE 

H=( 2 .+u.525*SQRT (REYNUM) ) *GSPH*GV I S/ ( 1 . 34*RAD ) 
Z=DDMDT*OVSPU( TEMBAR ) / ( H*DARE A ) 

IF ( Z-l • E— 4  )  2CU,3oO,3uU 
200  CONTINUE 
Z=1 • F-4 
GO  TO  5  no 
300  CONTINUE 

IF(25.-Z)  4  J  0 »  5  0  O’ » 5  0  0 
400  CONTINUE 
Z  =  2  5  o 

500  CONTINUE 

HDNET=H*uAREA*(GTEM-TFM)*Z/(EXP(Z)-l.  ) 

HO=HDNET -DDMDT * OHEVAP ( T CM) 

DV0L=4. /3.*3. 1415926*RAD**3 
HCAPD  =  ODENS(TEM) *DVOL*OLSPH<  TEM) 

DLTIM=DLPOS/VEL 
DLTEM= ( HD/HCAPD ) *DLT  IM 
I F ( TEM+DLTEM )  700,700,601 
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DLTEM=7. 

CONTINUE 
DLT IM=DLP0S/VEL 
C0DRAG  =  27./REYN’jy**  CJ.B4) 

ACCEL = ( 3./ 8« ) #COURAG* (GVlL-VtL )  *AoS < GVEL-VLL ) 

X  *(G*TM#GPRE/(UGCON*TEMtSAK)  )/(Ol>ENS<  TEM ) *R AO ) 

DLVEL=ACCEL*DLT!M 
IF  (VEL  +  DLVEL  )  7^.700,602 
?  CONTINUE 

DLRAD  =  -DDMDT*DLT  IM/  ( UAR EA#0LENS (  TEN*)  ) 

X  -RAD# ( ODENS ( TEM+DLT EM ) -ODENS ( TEM ) ) / ( 3 . #ODENS ( T EM ) ) 


I F ( RAD  +  DLR AD )  700, 70f 
603  CONTINUE 
60S  CONTINUE 

I  F ( ABS ( DLTEM ) .GE.10. 
GO  TO  750 
700  CONTINUF 
POP=0, 

RAD=0. 

DLRAD=0. 

T  E  M  =  0  . 


GO  TO  100 


DETEK=0. 

TIM=0. 

DLT  IM='/» 

X DO MOT  =  0. 

DDMDT  =  0 .  ' 

VEL=  0 • 

DLVEL=0. 

750  CONTINUE 
XPOS=POS 
XHEP=HEP 

xshint=shint 

XPOP=POP 
XRAD=RAD 
XTEM=JEM 
X'IM=T IM 
X  VFL  =  VEL 

FX=ODENS ITEM ) #  RAD##3  ■ 

RAD--,{AD+DLRAD 
TEM=TFM+DLTEM 
FY  =  ODENS  (  TE'-‘ )  *RAD*#3 
T IM=T IM+DLT IM 
P0S  =  P0S  +  DLPOS 
VEL=VEl+DLVEL 

SH I  NT  =  SH  I  NT  +  OLSPh  (  TEM  )  *l>L  TfcM 
DT=DLT IM 

HEP=HEP+ ( SHINT+Clil VAP( TEM) >#1.33  3  33*3. 1 4  1  59*AbS ( FX-FY ) #POP/DT 
1*1 ./FLOAT ( NSUBDV ) 

WBNUM  =  G'.','TM#GDPr  #  (  GVrL-VFL  )**2*?.*(PAD  ) 

X  /  (  UGCON#  T  E’-‘13  AR*OS  TF  N  (  TEM)  ) 

IF  O/FBRCR-WPN'JM)  8  00  ,  SCO  .  1 00"' 

800  CONTINUE 

RAD=RAD/FDSriNC#*.333333 

pop=pop*odshno 

0MEGA=RAD#«2/XRAD*«2*0MEGA 
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DDf-'DT  =  Of-T  C-A 
lOi.T  CONTINUE 
1C  )1  CONTINUE 

WRITE(6,2222)  NXTPT, NSUBDV 

2222  FORMAT  (1X2IB) 

ODFORI NXTPT  )  =  DVOL*UDENS <  T EM ) #DLVtL/ < DL T I M*GRCON ) 
ODFOVI NXTPT )=ODFOR< NXTPT )*VEL 
OOHES ( NXTPT ) =  HEP 
HEP  =  0 • 

ODHET( NXTPT )=HDNET 
ODPOP(NXTPT)=POP 
ODRAD(NXTPT) =  RA0 
ODTEM ( NXTPT ) =T  EM 
ODT I M( NXTPT ) =T IM 
ODVEL (NXTPT)=VEL 
RETURN 

99  IF(J.EQ.l)  TEM=TEM-ABS(DLTEM) 

100  NN=J*2-3 

DLPOS=DLPOS/2.0 

N5UBDV=NSUBDV*2 

IF(NN.LT.O)  GO  TO  113 

POS^XPOS 

POP=XPOP 

RAD=XRAD 

TEM=XTEM 

T  I M  =  X  T  !  M 

VEL  =  XVFL 
HEP=XHEP 
SHINT=XSHINT 
GO  TO  126 
113  NN= 1 

GO  TO  126 
END 


SlbFTC  COPROP  DECK 

SUBROUT I NE  CGPROP 

DIMENSION  CGDfc'M  I  1 29 ) , CGFLO ( 1 2 9 > .CGPOSl 129) ,CGPRE( 1 29 ) ,CGT EM ( 129 ) , 
1  CGVEL ( 1 29) ,CGWMO( 129) 

DIMENSION  FDRAT( 129) , r I  RAT  <  129) 

DIMENSION  FBFLO ( 129>,FBFOR< 129) ,FBFOV( 1 29  )  ,  FBFIES  (  129  ) *FBHET ( 129) » 

1  FBPOP  (  129)  , FBVEL (  129)  » 

2  FDFOR ( 129)* FDFOVt 129 ) .FUHtSl 129 ) ,huHtT ( 129) ,FDPOP( 129) , 

3  FDPOSf 129) i FDR AD ( 129 ) ,FDT EM ( 1 29 ) , P'DT I M ! 129 ) ,FDVEL(  129) » 

4  FMFLOt 129 ) ,FMFOR( 129 ) »FMFOV< 129 ) *  FMHES ( 129) .FMHETl  129 ) , 

6  .  FMPOP  (  129  ). 

DIMENSION  FDCONB ( 33 ) ,FDPOSB(33) .FDRADB ( 33 ) »  FDJEMB ( 33 ) ,FDVELb(33) 
DIMENSION  OBFLOf 129) ,OBFOR< 129) ,OBFOV( 129) ,OBHES( 129) *OBHET ( 129) * 

1  OBPOP(129),OBVEL(129), 

2  ODFOR  (  129  )  .ODFOVI  129  )  »Ol)HES(  129  )  »ODHET  (129)  ,ODPOP(  129)  , 

3  ODPOS  (  129  )  ,ODRAD(  129)  ,ODTEM( 129  )  .OUT  I M (  129)  ,GDVEL( 129)  . 

4  OMFLO  (  129)  *  OMFOR  (  1291  ,OMFOV(  129)  »OMHES(  129)  .OMHEK129)  ♦ 

5  OMPOP (129) 

DIMENSION  ODCONB ( 33 ) ,0DP0SB(33) .0DRADBI33) ,ODTEMB(33) ,ODVELB(33) 
DIMENSION  PAGET  I  (  12) .PROGTI ( 12) 

DIMENSION  SCGDFN (  129) ,SCGFLO(129) tSCGPREl 129) , SCOT  EM (129) , 

1  SCGVEL (  129  ) 

COMMON  CGDEN,CGFLO,CUPOS»CGPRE,CGT EM, CGVEL »CGWMO 
COMMON  EDR AT . £  I  RAT 

COMMON  FBFLO, FBFOR.FbFOV.FBHES.FBHET  » F BPOP » Fb VEL » 

1  FDFOR, FOFOV.F DUES, FUHET , FUPOP » FDPUd , FDR AD , FDTEM, FDT IM.FDVEL, 

2  FMFLO.FMFOR ,FMF0V, FMHES, FMHET, FMPOP 
COMMON  FDCONB, FDPOSB, FDRADB, FDTEMb.FDVE LB 
COMMON  OBFLO,OBFOR,OBFOV,OBHES,OBHET,OBPOP,OBVEL, 

1  ODFOR, ODFOV,ODHES,ODHET,ODPOP, ODPOS, ODR AD, ODTEM.ODT I M.ODVEL, 

2  OMFLO , OMFOR , OMFOV ,OMHES ,OMHET .OMPOP 
COMMON  ODCONB, ODPOSB .ODRADB ,ODTEMB , ODVELB 
COMMON  PAGET! .PROGTI 

COMMON  SCGDEN, SCO FLO.oCGPRE, SCOT EM, SCGVEL 

COMMON  chambl,conczl,dinozl,fdshno,fjfloi .FJVELI .fwtmol.genths, 

1  GWTMOL.ODSHNO.CJFLOI , OJVEL I , OWTMOL , V ARM Ax , WEBRCR 

COMMON  NOFORS ,NOI TRS .NOODRS .NOOPTS 
COMMON  LTIN.LTCUT 

EQUIVALENCE  ( NOFPTS.NUGPTS, NOOPT  S , NOPNT  S ) 

DO  230  1=3  .  NOPN  T  S 

IF ( FBFLO ( I ) -ToFLOl l-l ) )  21u,22u,21u 
210  CONTIN'JF 

EDRAT ( 1-1  )  =  ( OB  FLO ( I  )-ObFLO( 1-2) )/( FBFLO ( I ) -FbFLO ( 1-2)  ) 

GO  TO  230 
220  CONTINUF 

EDRAT (  1-1  )=0. 

IF ( FRF  LO (  I ) .EQ.FJFLOI )  LOR  AT { 1-1 ) =luU. 

IF (ORFLOt  I ) .EU.OJFLOI )  EDR  h I ( I  ) =u , 

230  CONTINUE 

FDR  AT ( 1  ) =0, 

EDRAT (NOPNTS)= EDRAT (NOPNTS-1 ) 

DO  260  1  =  1  .NOPNTS 
I  F  (  F  JFlO  !  rSrLO(D)  240,250,240 
240  CONTINUE 

F I R AT ( I , = (OJFLOI-OBFLOI I ) ) / ( F JF LO I -F dF LO ( 1 ) ) 


100 


go  to  2  6: 

2  5)  CONTINUE 

E  I  RAT ( I ) =130, 

IF(0JFL0I .EO.OBFlOI I ) )  EIRATl I )=, . 

260  CONTINUE 

COMMON/ROUORY/PA,TA ,v a,  ini t pt  »EA 
LLL= 1NITPT+1 

13  DO  200  I =LLL  *  NOPNT S 

CGFLOU )=AMAX1 (U. , ( FJFLOI -FBFLOf I ) +OJFLO I -OBFLO ( II ) ) 

S  C  G  F  L  0 ( I )=CGFLO< I ) 

CGVEL ( I ) =CGFLO ( 1 ) * 1 545 .*CGT EM ( I ) / ( CGWMOL ( CGTEM1 I ) , E I R AT ( I ) ) * 
1  CGPREI I )*CHAREA(FDPOS( I ) ) ) 

CGDENI I  1 =CGFLO ( I ) / ( CHAREA ( FDPOS ( I )  )*CGVEL( I  ) ) 

SCGDEN { I ) =CGDEN ( l ) 

200  CONTINUE 
KK=LLL+1 

DO  400  NXTPT=KK,NOPNTS 
DEN=CGDEN(NXTPT-1) 

POS  =  FDPOS (NXTPT-1 ) 

PRE  =  CGPRE (NXTPT-1 ) 

TEM=CGTEN (NXTPT-1 ) 

VEL  =  CGVEL (NXTPT-1 ) 

CALL  CGADV(NXTPT*DEN,POS,PRE.TEM.VEL) 

4u0  CONTINUE 

CALL  MAKW0N(CUD.DUM,N9  ) 

IF  (N9.LT.0)  GO  TO  13 
DO  500  1=1 .NOPNTS 

CGWMO ( I ) =CGWMOL ( CGTEM ( I ) » E I  RAT  (  I  ) ) 

500  CONTINUE 
RETURN 
END 


SI3FTC  CGAOV  DECK 

SUBROUTINE  CGADV< NXTPT . DEN .POS ,PRE , TEM, VEL ) 

DIMENSION  CGDENI 129) ,CGFLO( 129) ,CGPOS( 129) >CGPRE( 129 ) »CGT EM( 129 ) » 
1  CGVEL( 129) ,CGWMO( 129) 

DIMENSION  EDRAT  < 129) *EIRAT  C 129) 

DIMENSION  FBFLOI129) ,FBFOR( 129) ,FBFOV( 129) »FBHES(  129) »FBHET( 129) » 

1  FBPOP( 129) ,FBVEL( 129) ♦ 

2  FDFOR ( 1 29 ) » FDFOV ( 129 ) .FDHES ( 129 ) * FDHET ( 129 ) « FDPOP ( 1 29 ) » 

3  FDPOS< 129) ,FDRAD< 129 ) »FDT  EM( 129 ) , r DT IM( 129 ) »FDVEL (129), 

4  FMFLO ( 1 29 ) *  FMFORl 129 ) ,FMFOV( 129)  »FMHES( 129) »FMHET (129) . 

5  FMPOP (129) 

DIMENSION  FDCONB ( 33 ) .FDPOSB ( 33 ) »FDRADB( 33) , FDTEMB ( 33 ) .FDVELBI 33 ) 
DIMENSION  OBFLO(129) ,OBFOR( 129) ,OBFOV( 129) *OBHES(  129 ) ,OBHET( 129) . 

1  OB POP ( 129 ) * OB VEL ( 129 ) * 

2  ODFOR ( 1 29 ) »ODFOV( 129 ) *ODHES( 129) *ODHET( 129) »ODPOP(129) » 

3  ODPOS( 129 ) .ODRADl 129 ) »ODTEM( 129 ) *QDT IM( 129 ) .ODVEH 129 ) * 

4  OMFLO( 129) »OMFOR( 129) *OMFOV( 129)  .OMHESt 129) .OMHET (129) » 

5  OMPOP ( 129 ) 

DIMENSION  ODCONB( 33 ) *ODPOSB ( 33 ) *ODRADB( 33) .ODTEMB ( 33 ) ,ODVELB( 33 ) 
DIMENSION  PAGET  1(12) *PROGT 1(12) 

DIMENSION  SCGDEN( 129) *SCGFLO( 129 ) . SCGPRE ( 129 ) »SCGTEM ( 129 ) * 

1  SCGVEL ( 129 ) 

COMMON  CGDEN»CGFLO*CGPOS  *CGPRE  *CGTEM»CGVEL  *CGwMO 
COMMON  EDRAT, El  RAT 

COMMON  FBFLO,FBFOR,FBFOV,FBHES,FBHET,FBPOP,FBVEL» 

1  FDFOR, FDFOV, FDHES, FDHET, FDPOP, FDPOS, FDR AD, FDTEM, FDT I M.FDVEL, 

2  FMFLO, FMFOR.FMFOV.FMHES.FMHET, FMPOP 
COMMON  FDCONB, FDPOSB, FDRADB, FDT EMB,FDVELB 
COMMON  OBFLO,OBFOR,OBFOV,OBHES,OBHET,OBPOP,OBVEL» 

1  ODFOR  *ODFOV ,ODHES ,ODHET ,ODPQP ,ODPOS ,ODRAD  »ODTEM,ODT I M,ODVEL , 

2  OMFLO  »OMFOR ,OMFOV ,OMHES .OMHET , OMPOP 
COMMON  ODCON8 , ODPOSB , ODRADB .ODTEMB , ODVELB 
COMMON  PAGET  I , PkOGT I 

COMMON  SCGDEN.SCGFLO, SCGPRE. SCGTEM, SCGVEL 

COMMON  CHAMBL , CONOZL , D I NOZL , FDSHNO , F JFLOI , FJVEL I , FW f MOL .GENTHS, 

1  GWTMOL.ODSHNO.OJFLOI ,OJVEL I .OWTMOL , VARMAX , WEBRCR 

COMMON  NOFDRS.NOITRS, NOODRS ,NOOPT S 
COMMON  LT IN.LTOUT 

EQUIVALENCE  (NOFPTS.NOGPTS.NOOPTS.NOPNTS) 

COMMON/TEST/NA ,NB»PM IN ,PMAX 
COMMON/ MAK/ AM ACH( 129 ) 

COMMON /FORCE/N9 

COMMON/BOUDRY/PA.TA.VA, INI TPT  , EA 
NSUBD=8 
GRCON=32.2 
HETMF=77B. 

UGC0N=1545. 

SDEN=DEN 

SPOS=POS 

SPRE=PRE 

stem=tem 

SVEL=VEL 
KZ 1 =KZ 1 +1 
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IFlFLOX.LE.l.F-10  .AND.FDPOStNXTPT ) oLTeCHAMBL)  NSUBD^l 
IF(KZ1.LEoN0PNTS/30)  NSUBD  =6000/N0PNTS 
I F (  KZ1.GE.NOPNTS-INITPT-2-NOPNTS/30)  NSUBD=6000/NOPMTS 
IFlNXTPT .EQ.NOPNTS)  KZ1  =  0 
100  CONTINUE 
DEN=SDEN 
POS=SPOS 
PRE=SPRE 
TEM=STEM 

vel=svel 

DLP0S= ( FDPOS ( NX  TPT ) -POS ) /FLOAT ( NSUBD ) 

DO  1000  J=1«NSUBD 
IF(NTEST.EQ.l )GO  TO  13 
XDEN=DEN 
XPRF=PRE 
XEMACH=FMACH 
XTEM=TEM 
XVEL=VEL 

FRAC1= ( POS-FDPOS ( NXTPT-1 )) / ( FDPOS <  NXTPT ) -FDPOS( NXTPT-1 ) ) 

FRAC2=(  (  POS+DLPOS ) -FDPOS (NXTPT-1 ) ) /( FDPOS ( NXTPT )-FDPOS( NXTPT-1 > ) 
FLO=SCGFLO< NXTPT-1 )+FRACl*(SCGFLO(NXTPT)-SCGFLO( NXTPT-1 ) ) 

FLOX=< SCGFLO (NXTPT )-SCGFLO( NXTPT-1 > )/( FDPOS ( NXTPT 1-FDPOS ( NXTPT-1 ) > 
FFLOX  =  -( FBFLOI NXTPT )-FBFLO( NXTPT-1 ) ) / ( FDPOS ( NXTPT l-FDPOS ( NXTPT-1' ) ) 
OFLOX=-(O0FLO( NXTPT )-OBFLO< NXTPT-1 ))/(ODPOS( NXTPT )-ODPOS( NXTPT-1 > ) 
BFOR=( 1.-FRAC2 )*(FBFOR( NXTPT-1 I+OBFOR (NXTPT-1 ) )  + 

X  FRAC2*M  FBF0R  (  NXTPT  ) +OBFOR ( NXTPT  )  ) 

BFOV=( 1.-FRAC2 )*(FBFOV (NXTPT-1 )+OBFOV( NXTPT-1 ) )  + 

X  FRAC2*(FBF0V< NXTPT )+OBFOV(NXTPT  > ) 

BHES=FBHES( NXTPT )+OBHES( NXTPT) 

BH-T  =  ( 1 .-FRAC2 ) * ( FBHET ( NXTPT-1 ) +OBHET ( NXTPT-1 ) )  + 

X  FRAC2*( FBHET ( NXTPT ) +OBHET ( NXTPT ) ) 

BFVEL= ( 1 .-FRAC2 ) *FBVEL ( NXTPT-1 )  +FRAC2*FBVEL ( NXTPT ) 

BOVEL= ( 1 0-FRAC2 ) *OBVEL ( NXTPT-1 )  +FRAC2#0bVEL ( NXTPT ) 

AREA1=CHAREA ( POS ) 

AREA2=CHAREA ( POS+DLPOS ) 

FL0W1=FL0 

IF(  FLOV/1  )  110,150,110 
110  CONTINUE 

FL0W2=FL0+FL0X*DLP0S 

EIR1 =(1.-FRAC1)*EIRAT( NXTPT-1 )+FRACl*EI RAT (NXTPT) 

E I R2  =  ( 1. -FRAC2 ) *E IRAT ( NXTPT-1 )+FRAC2*E I  RAT (NXTPT) 
ENTS1=CGENTS(FLOX.FFLOX,OFLOX) 

D=DEN 

V=FLOW2/(AREA2»D) 

P=PRE-D»V*(V-VEL)/GRCON-(8FOR/AREA2)*DLPOS 
T  =  TEM 
HH»0. 

CAUL  DP ( E IR1 ,HH»TEM) 

XXl=FLOWl  *(HH+VEL«VEL/(2.*GRC0N*HETME) ) 

XX2=  (BHES-BHET+ENTS1«FL0X)  *DLP0S 

XX= (XX1+XX2 ) /FL0W2 
DO  400  K=1 ,50 
SAVED=D 
SAVEP=P 
SAVET=T 


SAVEV=V 

V=FL0W2/ ( A RE A 2* SAVE D ) 

G-ENTH  =  XX-SAVEV<>»2/ ( 2 . O^GRCON&HETME ) 

TT  =  0o 

CALL  DPI EIR2,GGENTH»TT  ) 

T  =  TT 

300  CONTINUE 

P=PRE-SAVED»SAVEV» I SAVEV-VEL ) /GRCON- f  BF0R/AREA2 ) *DLPOS 
D=SAVEP»CGWMOL  <  SAVET  » E l R2 ) / < SAVET»UGCON ) 

KQ=>K 

IF(KQ»EQo1  )  50T0400 
DEVD=ABS< (D-SAVED) /SAVED) 

DEVP=ABS( <P-$AVEP)/SAVE P) 

DEW  =  ABS<  (V-SAVEV)ZSAVEV) 

DEVMX  =  AMAX1 1 DEVD  * AMAX 1 1 DEVP  * AMAX1 1 DEVT  »DEVV ) )  ) 

I F ( 1.0E-3.GE.DEVMX)  GO  TO  500 
400  CONTINUE 
500  CONTINUE 

OLDM=EMACH 

IF(NXTPT.LT.NOPNTS/3)0LDM=0. 

EMACHaV/SQRT ( CGSPHT ( T ♦ E I R2 ) »GRCON« 1 545. * T / l CGSPHT ( T . E I R2 ) #CGWMOL ( 
1  T.EIR2J-1.958) ) 

DLMAK=FMACH-OLDM 
LOGICAL  POINT 

POlNTaFDPOS(NXTPT-l  ),GT.CHAMBL 

IFIEMACH. LE. 1.0. AND. < DLMAK.GE.O.. OR.  .NOT.  POINT) »G0  TO  900 
IF(NXTPT.GE.NOPNTS)  GO  TO  1025 
WRITE (6*1 )  OLDM,EMACH,NXTPT,J 

1  FORMAT (10X9HOLD  MACH=  F6 . 3 , 5X5HMACH=  F6. 3  * 5X6HNXTPT= 1 4 » 5H  S=I3) 

NA*1 

IF(NB.E0.1)G0T0135 

PA»PA»1.05 

GOTO? 

135  PMIN=PA 

PAo (PMIN+PMAX) /2,0 
7  CGPREI INITPT)=PA 

CGPREI IN1TPT+1 )=PA 

KZ1=0 

N9«-l 

EXTERNAL  CGPROP 
CALL  SNEAKY (CGPROP) 

651  CONTINUE 
900  CONTINUE 
DEN»D 
PRE=P 
TEMnT 
VFi  -V 

POS=POS+DLPOS 

C  IF  IT  IS  DESIRED  TO  INSERT  A  SUBSTEP  HALVING  OPTION  INTO  THIS 
C  SUBROUTINE,  JUST  INSERT  THE  TEST,  AND  IF  THE  SUBSTEP  LENGTH  IS 
C  CHANGED,  MAKE  A  SUBSEQUENT  TRANSFER  TO  100,  OTHERWISE  GO  ON  AS  USUAL, 
C  THIS  15  A  GOOD  PLACE  TO  INSERT  THE  TEST  IF  ONE  IS  NEEDED. 

1000  CONTINUE 

CGDEN ( NXTPT ) =DEN 
CGPRE(NXTPT) »PRE 
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CGTEM(NXTPT)=TEM 

H=GGENTH 

CGVEL(NXTPT)=VEL 
AMACH ( NXTPT ) =EMACH 
HO=H+V*«2/50100. 

T0=0, 

CALL  DP< EIR2*H0,T0) 

COMMON/TOHO/CGTO( 129 ) »CGHO( 129) 

CGTO(NXTPT)=TO 
CP=CGSPHT ( T  »  E I R2 ) 

GAMMA=CP/(CP-< le986/CGWMOL( T  »EIR2 ) ) ) 

PEX=GAMMA/ < GAMMA-1 . ) 

PO=P»<  TO/T)**PEX 
CGHO< NXTPT) =P0 

WRITE <6, 881 )NXTPT.P,V.T.FL0.EIP1 iGGENTH.EMACH 
6  FCRMAT (1X«3I4«5E16>8) 

RETURN 

1025  NTEST  =0 

IF(NXTPT.NE.NOPNTS)NTEST=l 

WRITE(6»6)NXTPT  *NSUB0 o J *XTEM »XPRE »XDEN»XVEL »XEMACH 
881  FORMAT  (lXI4t7E1.5»8i 
DO  5  N=NXTPT»NOPNTS 
CGDEN(N) =XDEN 
CGTEM' N ) =XTEM 
CGPRE ( N ) =XPRE 
CGVEL ( N ) =XVEL 
AMACH(N)=XEMACH 
5  CONTINUE 
RETURN 
150  D=DEN 
P  =  PRF 
V=VEL 
T  =  TEM 
GO  TO  900 
13  NTEST  =  0 
RETURN 
END 
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SIBFTC  MA<WON  DECK 

SUBROUTINE  KAKWOM VEL  ,  TEM.N9) 

DIMENSION  CGDEN ( 129) .CGFL0(129) , CGPOS! 129) .CGPRE! 129) »CGT EM(  129) » 
1  CGVEL  (129)  *  CG'a'MO  (129) 

DIMENSION  EDRATI 129) *  E  I  RAT { 129) 

DIMENSION  FbFLO ( 129 ) » FBFOR ( 1 29) ,FbFOV( 129) ,FdHES( 129 ) ,FbHET( 129) » 

1  .  FBPOPI 129) ,FBVEL( 129)  , 

2  F DFOk  (  129)  »  FDFOV  (.129  )  tFDHESI  129)  ,FUHET(  129)  *  FDPOP  1129)  * 

3  FDPOSI 129) ,FDRAD< 129)  ,  FDT EM < 1 29 ) , FDT I M ( 1 29 ) ,FDVEL(129) , 

A  FMFLOI 129) *  FMFOR  M  29 )  ,FMFOV< 129 ) , FMHES ( 129) ,FMHET ( 129) » 

5  FMPOP (129) 

DIMENSION  FDCONBl 33) ,FDPOSB( 33) *FDRADB( 33) <  FDTEMB ( 33 ) ,FDVELB(  33 ) 
DIMENSION  OBFLOI 129) ,OBFOR< 129) ,OBFOV( 129) ,OBHES( 129 ) ,OBHET ( 129 '  * 

1  OB POP ( 129) *OBVEL( 129)  , 

2  ODFORI 129) ,ODFOV( 129)  ,ODHES< 129) ,ODHET( 129)  .ODPOP1129) , 

3  ODPOSI  129)  ,  ODRAD  (  129  )  ,ODTEM( 129  )  »ODT  IM(  129)  ,ODVEL(129)  . 

A  OMFLO (129) * OMFOR ( 129)  ,OMFOV< 129) ,OMHES( 129) ,OMHET<129) * 

5  OMPOP (129)  . 

DIMENSION  ODCONBI33) .ODPOSBI33) .ODRADBI 33) »ODTEMB ( 33 ) *ODVELB ( 33 ) 
DIMENSION  PAGETI ( 12) »PROGT I ( 12) 

DIMENSION  SCGDEN! 129) ,SCGFLO( 129) ,SCGPRct 1 29) ,SCGTEM( 129) , 

1  SCGVEL I  129 ! 

COMMON  CGDEN, CGFLO. CGPOS. CGPRE.CGTEM. CGVEL  »CGWMO 
COMMON  EDRAT , E I  RAT 

COMMON  FBFLO, FBFOR, FBFOV , FBHES , FBHE  T , FBPOP , FBVEL, 

1  FDFOR , FDFOV , FDHtS, FDHET , FDPOP ,FDPOS , FDRAD , FDTLM, FDT I M , FDVEL , 

2  FMFLO, FMFOR , FMFOV, FMHES , FMHET , FMPOP 
COMMON  FDCONB, FDPOSB , FDR ADB, FDTEMB, FOVE LB 
COMMON  OBFLO,OBFOR,OBFOV,OBHES,ObHET,OBPOP,ObvEL» 

1  ODFOR,ODFOV,ODHES,ODHLT,ODPOP,ODPOS»ODRAD,ODTLM,ODTIM,ODVEL, 

2  OMFLO, OMFOR, OMFOV.OMHES.OMHET, OMPOP 
COMMON  ODCONB, ODPOSB ,ODRADB,ODTEMB,ODVELB 
COMMON  PAGETI .PROGTI 

COMMON  SCGDEN , SCGFLO , SCGPRE  ,  SCGT EM , SCGVEL 

COMMON  CHAMbL .CONOZL , D I N02L , FDSHNO , F JFLO I , FJVEL I , rWTMOL .GENTHS , 

1  GWTMOL.ODSHNO.OJFLCI .ojveli ,owtmol,varmax,webrcr 
COMMON  NOFDRS  »NOI TRS .NOODRS .NOOPTS 
COMMON  LTIN.LTOUT 

EQUIVALENCE  (NOFPTS.NOGPTS, NOOPTS. NOPNTS) 

COMMON/TEST/NA ,NB,PKIN,PMAX 
COMMON/BOUDRY/  PA , TA  , VA  ,  I N I TPT 
COMMON  /MAK/  AMACH (129) 

REAL  MACH 

MACH= AM ACHI NOPNTS) 

PAD=PA 

IF(MACH.LT,,97.0R.MACH.GT.1.03)  GO  TO  10 
N9  =  + 1 

9  WRITEI6.20)  MACH, PAD, PA 

20  FORMAT (10X,1PC)5.7,5X,E15.7,5X,E15.7) 

RETURN 

10  CONTINUE 

NB=  1 

IF(NA.EQ.l)  GOTO  135 

PA=MACH»PA 

GO  TO  7 
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P‘.'AX  =  PA 

P  A  =  (  PMIN  +  P'-’AX  )  /  2 . 0 
CGPPr (  !  N  I  TPT  t  =PA 
CGPRE ( INITPT+1 ) =PA 
N9  =  -  1 
GO  TO  9 
END 


? ;bftc  op 

SUBROUTINE  DP(60F,H,T) 

REAL  M 

COMMON/GOFTH/  TAI3LEI8.39I 

DIMENSION  SLOPE (2)  ,CONST<2) 

DO  10  <=2,32765 

1J  IF (TABLE ( 1 »K ) .GE.GOE )  GO  TO  20 

20  K=K- 1 

9  FRACTMGOF-TABLt  ( 1,0  )/ (TABLE!  1  ,<  +  l  i -TABLE  (1,0  ) 

DO  50  J  =  1  ,2 
JJ=K+ J-l 

SLOPE ( J)  =  ( TABLE ( 4,  JJ) -TABLE!  2,  JJ) ) / ( T ABLE ( 5 , JJ ) -T ABLE (  3 , J J ) ) 
50  CONST ( J ) =-  (  (SLOPE! J)*T ABLE (3, JJ) ) -TABLE ( 2  ,  JJ ) ) 

B=FRACT*( CONST ( 2  )  -CONST ( 1 ) l+CONST ( 1 ) 

M  =  FRACT* (SLOPE ( 2  ) -SLOPE! 1 ) J+SLOPE ( 1 ) 

I F ( Ho  EQ . 0 • )H= ( T-Bl/M 
IF(T„EQ.O.)T=M*H+B 
8  CONTINUE 

RETURN 
END 


S.IBFTC  TERP  DECK 

FUNCTION  TERP(XX,X,Y,NOXS) 

DIMENSION  X ( 2 1  , Y ( 2  ) 

IF(NOXS-I)  10,10,20 
10  CONTINUE 
TERP=Y ( 1 ) 

RETURN 
20  CONTINUE 

IF(XX-X(1))  100, 10U.200 
100  CONTINUE 

TERP=Y(l)  +  (  (Y(2)-Y(l) )/(X(2)-X(l)  >  ) * ( XX-X (  1)  ) 

RETURN 

200  CONTINUE 

I F ( X ( NCX S ) -X X  )  300 , 30  j, 400 

300  CONTINUE 

TERP  =  Y ( N0X5 )  +  (  (  Y  (  NOXS ) -Y ( NOXS-1 )  )  /  (  X  (  NOXS ) -X ( NOXS-1 )  )  ) 

X  * (XX-X (NOXS  )  ) 

RETURN 

400  CONTINUE 

DO  600  1=1 *NOXS 
I  F  (  X  (  I  )  -XX  )  500, 5 '0,6.' 0 
5C0  CONTINUE 
T  I  =  I 

600  CONTINUE 

TERP  =  Y( I  I  )+(  (Y(  I  I  +  1 )  -  Y ( II ) )/(X(  I  1  +  1) —X (II))  ) * ( XX-X ( II)) 

return 

END 
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S1BFTC  FHEVAP  DECK 

FUNCTION  FHtVAf'(TLM) 
DIMENSION  T ( 1 4 )  »  H ( 1 4 ) 

T ( 1 ) =540o 
T  (  2  )  =  6  6  0  , 

T ( 3 ) =71 0„ 

T(4) =760. 

T ( 5 ) =810. 

T(6) =860. 

T ( 7 ) =910. 

T<8)=960. 

T ( 9 1 =1010. 

T  C 10 ) =1060 • 

T(ll)=1110e 

T  ( 1 2 )  =1 160,"  — 

T(13)=1210» 

T(14)=1235. 

H  < 1 ) =150. 

H  (  2  )  =  142~« 

H  <  3 ) =137. 

H(4)=131. 

H ( 5 ) =125. 

H  C  6 ) =120. 

Hf 7) =11*. 

HIS) =106.  • 

HI  9) =98. 

HI  10 ) =88 . 

Hill ) =75. 

HI  12 ) =55. 

H 1 13 ) =25 . 

HI 14)=0. 

FHEVAP  =  TERPI1  EM,T,H,14) 

RETURN 

END 


$  I BFTC  OVPRE  DECK 

FUNCTION  OVPRE(TEM) 

0VPRE=144.*EXP I  11. 9584- 1476. 4912 /ITEM-3. 568)  ) 

RETURN 

END 


S1BFTC  CGEtiTS  DECK 

FUNCTION  CGENTSI  F.LOX  *FFLOX  *OFLOX  ) 
CGENTS=1.8*(-362.97*FFLOX-96.25#OFLOX) /FLOX 
•  RETURN 
END 
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PROPORTIONAL  OXIDIZER  EVAPORATION 

FROZEN  COMPOSITION 
(DATA  OF  REF.  4) 

-  r0=  10- 4  ft. 

- rQ  =  2  x  I0-4  ft. 

- * - r0=  3x  IO-4ft. 

-  - - r  =  4x  IO-4ft. 


\U 

I 


/  /  /T 

\!/Y 

W 


/ 


/ 


AXIAL  DISTANCE  (ftj 

FIG. 2.3  FUEL  DROP  TEMPERATURE  AND  VELOCITY 
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IG.  2.4  COMBUSTION  GAS  PROPERTIES 
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FIG.  2.5  DROP  TEMPERATURE  vs.  DISTANCE 
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FIG. 2.8  COMBUSTION  GAS  PROPERTIES  vs  DISTANCE 
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FIGURE  3.1  EXPERIMENTAL  ROCKET  ENGINE 


FIGURE  3.2  LIQUID  PROPELLANT  INJECTOR  HEAD 


FIGURE  3.3  500  LBF  INJECTOR  PLATE  BLANK 
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FIGURE  4.2  INJECTOR  SUMMARY:  GEOMETRY 
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FIGURE  4.3  AXIAL  PRESSURE  VARIATION  AND 
GRADIENT 
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FIGURE  4.4  AXIAL  PRESSURE  VARIATION 


FIGURE  4.5  AXIAL  PRESSURE  VARIATION 
AND  GRADIENT 


FIGURE  4.9  IMPINGING  INJECTOR  FOR  WAVE  PROPAGATION  EXPERIMENTS 
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FIGURE  4.14  WAVE  SLOPE  vs.  PRESSURE  GRADIENT 
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FIGURE  4.15  NORMALIZED  WAVE  SLOPE  vs.  PRESSURE  GRADIENT 
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